#include "wavefun_coeff.h"
#include <cmath>
using std::cout;
using std::endl;
using std::vector;
// ~~~~~~~~~some functions~~~~~~~~~~~~~~~//
int Wave_Fun::phase(int n) { return n % 2 ? -1 : 1; }

int Wave_Fun::delta(int a, int b) { return (a == b) ? 1 : 0; }

double Wave_Fun::ho_k(int n, int l, double b, double k) {
  // Becareful the phase = (-i)^(2n+l); Here only contain (-i)^(2n)
  double z = b * k;
  double temp =
      phase(n) * exp(-0.5 * z * z) * pow(z, l) * sqrt(2.0) *
      exp(0.5 * (gsl_sf_lngamma(n + 1) - gsl_sf_lngamma(n + l + 1.5))) *
      gsl_sf_laguerre_n(n, l + 0.5, z * z) * pow(b, 1.5);
  return temp;
}

double Wave_Fun::ho_r(int n, int l, double b, double r) {
  double z = r / b;
  double temp =
      exp(-0.5 * z * z) * pow(z, l) * sqrt(2.0) *
      exp(0.5 * (gsl_sf_lngamma(n + 1) - gsl_sf_lngamma(n + l + 1.5))) *
      gsl_sf_laguerre_n(n, l + 0.5, z * z) / pow(b, 1.5);
  return temp;
}

void Wave_Fun::gauleg(const double x1, const double x2, vector<double> &x,
                      vector<double> &w) {
  const double EPS = 3.0e-14;
  double z1, z, xm, xl, pp, p3, p2, p1;
  int n = x.size();
  int m = (n + 1) / 2;
  xm = 0.5 * (x2 + x1);
  xl = 0.5 * (x2 - x1);
  for (int i = 1; i <= m; i++) {
    z1 = 0.0;
    z = cos(Pi * (i - 0.25) / (n + 0.5));
    do {
      p1 = 1.0;
      p2 = 0.0;
      for (int j = 1; j <= n; j++) {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / (1.0 * (j));
      }
      pp = n * (z * p1 - p2) / (z * z - 1.0);
      z1 = z;
      z = z - p1 / pp;
    } while (std::abs(z - z1) > EPS);
    x[i - 1] = xm - xl * z;
    x[n - i] = xm + xl * z;
    w[i - 1] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
    w[n - i] = w[i - 1];
  }
}

void Wave_Fun::gauleg(const long double x1, const long double x2,
                      vector<long double> &x, vector<long double> &w) {
  const double EPS = 3.0e-14;
  long double z1, z, xm, xl, pp, p3, p2, p1;
  int n = x.size();
  int m = (n + 1) / 2;
  xm = 0.5 * (x2 + x1);
  xl = 0.5 * (x2 - x1);
  for (int i = 1; i <= m; i++) {
    z1 = 0.0;
    z = cos(Pi * (i - 0.25) / (n + 0.5));
    do {
      p1 = 1.0;
      p2 = 0.0;
      for (int j = 1; j <= n; j++) {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / (1.0 * (j));
      }
      pp = n * (z * p1 - p2) / (z * z - 1.0);
      z1 = z;
      z = z - p1 / pp;
    } while (std::abs(z - z1) > EPS);
    x[i - 1] = xm - xl * z;
    x[n - i] = xm + xl * z;
    w[i - 1] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
    w[n - i] = w[i - 1];
  }
}

void Wave_Fun::init(int n_max, int l_max, int lbarS, int EJC_max,
                    int TPE_RANK) {
  // if check_mpi == 1 check_mpi glbarS2_mpi PFKglbarS_mpi
  //
  if (TPE_RANK != WaveCom.rank()) {
    cout << " !!!! Tpe_rank != WaveCom.rank()  !!! " << endl;
    exit(0);
  }
  this->lbar_max = lbarS;
  double alpha_b = sqrt((1.0 * hbarc * hbarc) / (P_massave * hbar_omega));
  vector<double> k_t, w_t;

  // cout<<alpha_b<<endl;
  k_t.resize(k_mesh);
  w_t.resize(k_mesh);
  // this->gauleg(0.0,25.0,k_t,w_t);

  double *k_l;
  double *w_l;
  k_l = new double[k_mesh];
  w_l = new double[k_mesh];

  this->legendre_compute_glr(k_mesh, k_l, w_l);
  if (WaveCom.rank() == 0) {
    cout << "\t k_scale = " << this->k_scale << endl;
  }
  this->rescale(0.0, this->k_scale, k_mesh, k_l, w_l);
  for (int i = 0; i < k_mesh; i++) {
    k_t[i] = k_l[i];
    w_t[i] = w_l[i];
  }

  if (WaveCom.rank() == 0) {
    std::ofstream file_kmesh;
    file_kmesh.open("data/k_mesh.dat");
    for (int i = 0; i < k_mesh; i++) {
      file_kmesh << k_t[i] << "\t" << w_t[i] << endl;
    }
    file_kmesh.close();
  }
  // initial klam_state
  int klamb_max = lbarS + 4; // k^( (lbar + 4) + 1)
  klamb_struct klamb_state_t;
  klamb_vec.resize(lbarS + 1);
  for (int l = 0; l <= lbarS; l++) {
    klamb_vec[l].resize(klamb_max + 1);
    for (int A = 0; A <= klamb_max; A++) {
      klamb_vec[l][A].resize(klamb_max + 1);
      for (int B = 0; B <= klamb_max; B++) {
        klamb_vec[l][A][B].resize(klamb_max + 1);
        for (int C = 0; C <= klamb_max; C++) {
          klamb_vec[l][A][B][C].resize(klamb_max + 1);
          for (int D = 0; D <= klamb_max; D++) {
            klamb_vec[l][A][B][C][D] = -1;
          }
        }
      }
    }
  }
  klamb_state.resize(lbarS + 1);
  for (int l = 0; l <= lbarS; l++) {
    int index = 0;
    for (int A = 0; A <= klamb_max; A++) {
      if (A > l + 2)
        continue;
      for (int B = 0; B <= klamb_max; B++) {
        if (B > l + 2)
          continue;
        if ((A + B) < l || (A + B) > l + 2)
          continue;
        for (int C = 0; C <= klamb_max; C++) {
          if (C > (l + 4))
            continue;
          for (int D = 0; D <= klamb_max; D++) {
            if (D > (l + 4))
              continue;
            if ((C + D) < l || (C + D) > (l + 4))
              continue;
            if ((A + B + C + D) > (2 * l + 4))
              continue;
            klamb_vec[l][A][B][C][D] = index;
            klamb_state_t.lamb12_p = A;
            klamb_state_t.lamb12 = B;
            klamb_state_t.lamb3_p = C;
            klamb_state_t.lamb3 = D;
            klamb_state[l].push_back(klamb_state_t);
            index++;
          }
        }
      }
    }
  }
  if (WaveCom.rank() == 0) {
    int sum = 0;
    for (int l = 0; l <= lbarS; l++) {
      cout << " l =  " << l << "\tklamb_size = " << klamb_state[l].size()
           << endl;
      sum += klamb_state[l].size();
    }
    cout << "\tklamb_size sum = " << sum << endl;
  }
  //
  klamb_max = lbarS + 5; // k^(lbar + 5)
  Pwave_w.resize(k_mesh);
  Pwave_x.resize(klamb_max);

  for (int i = 0; i < klamb_max; i++) {
    Pwave_x[i].resize(k_mesh);
  }
  for (int i = 0; i < k_mesh; i++) {
    for (int j = 0; j < klamb_max; j++) {
      int k = j + 1;
      this->Pwave_x[j][i] = poww(k_t[i], k);
    }
    this->Pwave_w[i] = w_t[i];
  }

  Pwave_vec.resize(n_max + 1);
  for (int n = 0; n <= n_max; n++) {
    Pwave_vec[n].resize(l_max + 1);
    for (int l = 0; l <= l_max; l++) {
      if ((2 * n + l) > EJC_max)
        continue;
      Pwave_vec[n][l].resize(k_mesh);
      for (int i = 0; i < k_mesh; i++) {
        Pwave_vec[n][l][i] = this->ho_k(n, l, alpha_b, k_t[i]) * k_t[i];
      }
    }
  }

  Pwave_vecwx.resize(klamb_max);
  for (int m = 0; m < klamb_max; m++) {
    Pwave_vecwx[m].resize(n_max + 1);
    for (int n = 0; n <= n_max; n++) {
      Pwave_vecwx[m][n].resize(l_max + 1);
      for (int l = 0; l <= l_max; l++) {
        if ((2 * n + l) > EJC_max)
          continue;
        Pwave_vecwx[m][n][l].resize(k_mesh);
        for (int i = 0; i < k_mesh; i++) {
          Pwave_vecwx[m][n][l][i] =
              Pwave_x[m][i] * Pwave_vec[n][l][i] * Pwave_w[i];
        }
      }
    }
  }
  // calculate regulator vector
  this->regulator_cal();
  // calculate glbar vector
  // this->glbar_init(l_max);

  // initial kpk_state
  kpk_struct kpk_t;
  kpk_vec.resize(k_mesh);
  int index_kpk = 0;
  for (int i = 0; i < k_mesh; i++) {
    kpk_vec[i].resize(k_mesh);
    for (int j = 0; j < k_mesh; j++) {
      kpk_t.n_k = i;
      kpk_t.n_kp = j;
      kpk_state.push_back(kpk_t);
      kpk_vec[i][j] = index_kpk;
      index_kpk++;
    }
  }

  // kpk_struct kpk_t;
  // for (int i = 0; i < k_mesh; i++) {
  //   for (int j = 0; j < k_mesh; j++) {
  //     kpk_t.n_k = i;
  //     kpk_t.n_kp = j;
  //     kpk_state.push_back(kpk_t);
  //   }
  // }
  //// initial k123_state & k123_vec
  k123_struct k123_t;
  k123_vec.resize(k_mesh);
  int index = 0;
  for (int i = 0; i < k_mesh; i++) {
    k123_vec[i].resize(k_mesh);
    for (int j = 0; j < k_mesh; j++) {
      k123_t.n_k12 = i;
      k123_t.n_k3 = j;
      k123_state.push_back(k123_t);
      k123_vec[i][j] = index;
      index++;
    }
  }
  size_k123 = k123_state.size();

  /*if (WaveCom.rank() == 0 && evaluate_flag < 0) {
    this->PFKglbar0_build();
    //
    this->PFKglbar2_build();
    //
    // this->glbar_vec.clear();
    cout << "\t k_mesh = " << k_mesh << "\t p_mesh = " << p_mesh << endl;
    cout << "\t n_max = " << n_max << "\t l_max = " << l_max
         << "\t lbarS = " << lbarS << endl;
    cout << "nl_state size : " << nl_state.size() << endl;
    // cout<<"lbarS = "<<lbarS<<endl;
  } */

  /*
    Time_Count TC_2;
    TC_2.start();
    int lmax_legendrePlW = lbarS;
    int lbarS_G = lbarS;
    int lbarS_gg = lbarS_G;
    this->glbarS_init(lbarS_gg, lbarS_gg, lbarS_gg, lmax_legendrePlW);
    if (WaveCom.rank() == 0) {
      cout << "\t lbar_glbarS_max = " << this->lbar_glbarS_max << endl;
    }
    int lbarS_pp = lbarS_P;
    if (WaveCom.rank() == 0) {
      cout << "\t ?? \t good ??" << endl;
    }

    int mpi_run_flag = 1;
    vector<int> mpi_run_flag_gather(WaveCom.size());
    all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);

    vector<int> mpi_run_flag_gather(WaveCom.size());
    int mpi_run_flag = 1;
    if (evaluate_flag < 0 && MPI_TYPE == 1) {
      this->glbarS_vec_build_MPI();
    }
    if (evaluate_flag < 0 && MPI_TYPE == 2) {
      this->glbarS_vec_build_MPI_II();
    }
    if (WaveCom.rank() == 0) {
      cout << "glbarS_vec_build_MPI initial&calculate Time cost : " <<
    TC_2.sec()
           << "s" << endl;
      // exit(0);
    }/* - */

  init_flag = 1;
}

void Wave_Fun::init_single(int n_max, int l_max) {

  cout << "\t Wave_Fun::init_single  Start " << endl;
  evaluate_flag = -1;

  double alpha_b = sqrt((1.0 * hbarc * hbarc) / (P_massave * hbar_omega));
  vector<double> k_t, w_t;

  // cout<<alpha_b<<endl;
  k_t.resize(k_mesh);
  w_t.resize(k_mesh);
  // this->gauleg(0.0,25.0,k_t,w_t);

  double *k_l;
  double *w_l;
  k_l = new double[k_mesh];
  w_l = new double[k_mesh];

  this->legendre_compute_glr(k_mesh, k_l, w_l);
  if (WaveCom.rank() == 0) {
    cout << "\t k_scale = " << this->k_scale << endl;
  }
  this->rescale(0.0, this->k_scale, k_mesh, k_l, w_l);
  for (int i = 0; i < k_mesh; i++) {
    k_t[i] = k_l[i];
    w_t[i] = w_l[i];
  }

  if (WaveCom.rank() == 0) {
    std::ofstream file_kmesh;
    file_kmesh.open("data/k_mesh.dat");
    for (int i = 0; i < k_mesh; i++) {
      file_kmesh << k_t[i] << "\t" << w_t[i] << endl;
    }
    file_kmesh.close();
  }

  Pwave_w.resize(k_mesh);
  int klamb_max = 4;
  Pwave_x.resize(klamb_max);

  for (int i = 0; i < klamb_max; i++) {
    Pwave_x[i].resize(k_mesh);
  }
  for (int i = 0; i < k_mesh; i++) {
    for (int j = 0; j < klamb_max; j++) {
      int k = j + 1;
      this->Pwave_x[j][i] = poww(k_t[i], k);
    }
    this->Pwave_w[i] = w_t[i];
  }

  Pwave_vec.resize(n_max + 1);
  for (int n = 0; n <= n_max; n++) {
    Pwave_vec[n].resize(l_max + 1);
    for (int l = 0; l <= l_max; l++) {
      Pwave_vec[n][l].resize(k_mesh);
      for (int i = 0; i < k_mesh; i++) {
        Pwave_vec[n][l][i] = this->ho_k(n, l, alpha_b, k_t[i]) * k_t[i];
      }
    }
  }

  Pwave_vecwx.resize(klamb_max);
  for (int m = 0; m < klamb_max; m++) {
    Pwave_vecwx[m].resize(n_max + 1);
    for (int n = 0; n <= n_max; n++) {
      Pwave_vecwx[m][n].resize(l_max + 1);
      for (int l = 0; l <= l_max; l++) {
        Pwave_vecwx[m][n][l].resize(k_mesh);
        for (int i = 0; i < k_mesh; i++) {
          Pwave_vecwx[m][n][l][i] =
              Pwave_x[m][i] * Pwave_vec[n][l][i] * Pwave_w[i];
        }
      }
    }
  }
  // calculate regulator vector
  this->regulator_cal();
  // calculate glbar vector
  this->glbar_init(l_max);

  // initial nl_vec and nl_state for PFKglbar_vec
  //
  // initial kpk_state
  kpk_struct kpk_t;
  kpk_vec.resize(k_mesh);
  int index_kpk = 0;
  for (int i = 0; i < k_mesh; i++) {
    kpk_vec[i].resize(k_mesh);
    for (int j = 0; j < k_mesh; j++) {
      kpk_t.n_k = i;
      kpk_t.n_kp = j;
      kpk_state.push_back(kpk_t);
      kpk_vec[i][j] = index_kpk;
      index_kpk++;
    }
  }
  //// initial k123_state & k123_vec
  k123_struct k123_t;
  k123_vec.resize(k_mesh);
  int index = 0;
  for (int i = 0; i < k_mesh; i++) {
    k123_vec[i].resize(k_mesh);
    for (int j = 0; j < k_mesh; j++) {
      k123_t.n_k12 = i;
      k123_t.n_k3 = j;
      k123_state.push_back(k123_t);
      k123_vec[i][j] = index;
      index++;
    }
  }
  size_k123 = k123_state.size();

  // initial nl_vec and nl_state for PFKglbar_vec
  nl_struct nl_state_t;
  index = 0;
  nl_vec.resize(n_max + 1);
  for (int n12 = 0; n12 <= n_max; n12++) {
    nl_vec[n12].resize(n_max + 1);
    for (int n3 = 0; n3 <= n_max; n3++) {
      nl_vec[n12][n3].resize(l_max + 1);
      for (int l12 = 0; l12 <= l_max; l12++) {
        nl_vec[n12][n3][l12].resize(l_max + 1);
        for (int l3 = 0; l3 <= l_max; l3++) {
          nl_vec[n12][n3][l12][l3] = index;
          nl_state_t.n12 = n12;
          nl_state_t.n3 = n3;
          nl_state_t.l12 = l12;
          nl_state_t.l3 = l3;
          nl_state.push_back(nl_state_t);
          index++;
        }
      }
    }
  }
  this->PFKglbar0_build();
  //
  this->PFKglbar2_build();
  //
  if (WaveCom.rank() == 0 && evaluate_flag < 0) {
    // this->glbar_vec.clear();
    cout << "\t k_mesh = " << k_mesh << "\t p_mesh = " << p_mesh << endl;
    cout << "\t n_max = " << n_max << "\t l_max = " << l_max << endl;
  }
  init_flag = 1;
}

void Wave_Fun::regulator_cal() {

  // double hbarc_4 = pow(hbarc, 4);
  // double Cut_Lambda_4 = pow(Cut_Lambda, 4);
  double hbarc_2 = pow(hbarc, 2);
  double Cut_Lambda_2 = pow(Cut_Lambda, 2);

  F_reg_mat.resize(k_mesh, k_mesh);
  F_reg_vec.resize(k_mesh);
  for (int i = 0; i < k_mesh; i++) {
    double k12_2 = this->Pwave_x[1][i];
    F_reg_vec[i].resize(k_mesh);
    for (int j = 0; j < k_mesh; j++) {
      double k3_2 = this->Pwave_x[1][j];
      double x = (k12_2 + k3_2) * hbarc_2 / (2.0 * Cut_Lambda_2);
      x = -1 * pow(x, regulator_pow);
      F_reg_vec[i][j] = exp(x);
      // F_reg_vec[i][j] = 1.0;
      F_reg_mat(i, j) = F_reg_vec[i][j];
    }
  }
}

void Wave_Fun::glbar_init(int lmax) {

  // cal glbar_vec
  legendreP_x.resize(p_mesh);
  legendreP_w.resize(p_mesh);
  this->gauleg(-1.0, 1.0, legendreP_x, legendreP_w);

  legendre_Pl.resize(lmax + 2);
  for (int l = 0; l <= lmax + 1; l++) {
    legendre_Pl[l].resize(p_mesh);
    for (int u = 0; u < p_mesh; u++) {
      legendre_Pl[l][u] = gsl_sf_legendre_Pl(l, legendreP_x[u]);
    }
  }

  glbar_vec.resize(2);
  for (int i = 0; i < 2; i++)
    glbar_vec[i].resize(lmax + 2);

  // for glbar0
  for (int l = 0; l <= lmax + 1; l++) {
    glbar_vec[0][l].resize(k_mesh);
    for (int k3 = 0; k3 < k_mesh; k3++) {
      glbar_vec[0][l][k3].resize(k_mesh);
      for (int k3p = 0; k3p < k_mesh; k3p++) {
        double k3_x = this->Pwave_x[0][k3];
        double k3p_x = this->Pwave_x[0][k3p];
        glbar_vec[0][l][k3][k3p] = glbar_cal(k3_x, k3p_x, l, 0);
      }
    }
  }
  // for glbar2
  for (int l = 0; l <= lmax + 1; l++) {
    glbar_vec[1][l].resize(k_mesh);
    for (int k3 = 0; k3 < k_mesh; k3++) {
      glbar_vec[1][l][k3].resize(k_mesh);
      for (int k3p = 0; k3p < k_mesh; k3p++) {
        double k3_x = this->Pwave_x[0][k3];
        double k3p_x = this->Pwave_x[0][k3p];
        glbar_vec[1][l][k3][k3p] = glbar_cal(k3_x, k3p_x, l, 2);
      }
    }
  }
}

double Wave_Fun::glbar_cal(double k3, double k3p, int lbar, int g_index) {
  // using gsl_sf_legendre_Pl to calculate legendre polynomials

  double result = 0.0;
  // for glbar0
  if (g_index == 0) {
    double part1 = 0.0;
    for (int u = 0; u < p_mesh; u++) {
      double fx_t = k3p * k3p + k3 * k3 - 2 * k3 * k3p * legendreP_x[u];
      double fx = fx_t / (2.0 / 3.0 * fx_t + Mass_pi * Mass_pi);
      part1 += fx * legendre_Pl[lbar][u] * legendreP_w[u];
    }
    result = part1;
  } else if (g_index == 2) {
    double part1 = 0.0;
    for (int u = 0; u < p_mesh; u++) {
      double fx_t = k3p * k3p + k3 * k3 - 2 * k3 * k3p * legendreP_x[u];
      double fx = 1.0 / (2.0 / 3.0 * fx_t + Mass_pi * Mass_pi);
      part1 += fx * legendre_Pl[lbar][u] * legendreP_w[u];
    }
    result = part1;
  } else {
    std::cout << "wrong with g_index in wave_coeff.cpp" << endl;
  }
  return result;
}

void Wave_Fun::PFKglbar0_build() {
  int num = nl_state.size();
  PFKglbar0_vec.resize(num);

  cout << "Begin to initial PFKglbar0 " << endl;

  Pgs_Bar pgs;
  pgs.init(num, 30);
  for (int i = 0; i < num; i++) {

    pgs.PrintBar(i);

    PFKglbar0_vec[i].resize(num);
    int n12p = nl_state[i].n12;
    int n3p = nl_state[i].n3;
    int l3p = nl_state[i].l3;
#pragma omp parallel for schedule(dynamic, 1)
    for (int j = 0; j < num; j++) {
      // cout<<i<<"\t"<<j<<endl;
      int n12 = nl_state[j].n12;
      int n3 = nl_state[j].n3;
      int l3 = nl_state[j].l3;
      PFKglbar0_vec[i][j] = PFKglbar0_cal2(n12p, n12, n3p, n3, l3p, l3);
    }
  }
  pgs.EndBar();
  cout << " PFKglbar0 initial finished " << endl;
}

void Wave_Fun::PFKglbar2_build() {
  int num = nl_state.size();
  PFKglbar2_vec.resize(num);

  cout << "Begin to initial PFKglbar2 " << endl;

  Pgs_Bar pgs;
  pgs.init(num, 30);
  for (int i = 0; i < num; i++) {
    // waiting bar
    pgs.PrintBar(i);
    // calculate
    PFKglbar2_vec[i].resize(num);
    int n12p = nl_state[i].n12;
    int n3p = nl_state[i].n3;
    int l3p = nl_state[i].l3;
#pragma omp parallel for schedule(dynamic, 1)
    for (int j = 0; j < num; j++) {
      int n12 = nl_state[j].n12;
      int n3 = nl_state[j].n3;
      int l3 = nl_state[j].l3;
      PFKglbar2_vec[i][j].resize(3);
      for (int lambda1 = 0; lambda1 <= 2; lambda1++) {
        // cout<<i<<"\t"<<j<<"\t"<<lambda1<<endl;
        if (lambda1 == 0) {
          PFKglbar2_vec[i][j][lambda1].resize(1);
          int lbar = l3p;
          PFKglbar2_vec[i][j][lambda1][0] =
              this->PFKglbar2_cal2(n12p, n12, n3p, n3, l3p, l3, 0, lbar);
        }
        if (lambda1 == 2) {
          PFKglbar2_vec[i][j][lambda1].resize(1);
          int lbar = l3;
          PFKglbar2_vec[i][j][lambda1][0] =
              this->PFKglbar2_cal2(n12p, n12, n3p, n3, l3p, l3, 2, lbar);
        }
        if (lambda1 == 1) {
          int l_min = std::max(std::abs(l3 - 1), std::abs(l3p - 1));
          int l_max = std::min(l3 + 1, l3p + 1);
          if (l_min <= l_max) {
            int lbar_size = l_max - l_min + 1;
            PFKglbar2_vec[i][j][lambda1].resize(l_max + 1);
            for (int lbar = 0; lbar <= l_min; lbar++)
              PFKglbar2_vec[i][j][lambda1][lbar] = 0.0;
            for (int lbar = l_min; lbar <= l_max; lbar++)
              PFKglbar2_vec[i][j][lambda1][lbar] =
                  this->PFKglbar2_cal2(n12p, n12, n3p, n3, l3p, l3, 1, lbar);
          }
        }
      }
    }
  }
  pgs.EndBar();
  cout << " PFKglbar2 initial finished " << endl;
}

double Wave_Fun::PFKglbar0_cal2(int &n12p, int &n12, int &n3p, int &n3,
                                int &l3p, int &l3) {
  double result = 0.0;
  double temp_1 = 0.0, temp_2 = 0.0;
  vector<double> PF12_p(k_mesh), PF12(k_mesh);
  vector<vector<double>> PglbarP(k_mesh);
  double F_reg = 0.0;

  for (int k12 = 0; k12 < k_mesh; k12++) {
    temp_1 = 0.0;
    temp_2 = 0.0;
    for (int k3 = 0; k3 < k_mesh; k3++) {
      F_reg = F_reg_vec[k12][k3];
      temp_1 += Pwave_vecwx[0][n12p][0][k3] * F_reg;
      temp_2 += Pwave_vecwx[0][n12][0][k3] * F_reg;
    }
    PF12_p[k12] = temp_1;
    PF12[k12] = temp_2;
  }

  double P3p, P3;
  for (int k3p = 0; k3p < k_mesh; k3p++) {
    P3p = Pwave_vecwx[0][n3p][l3p][k3p];
    PglbarP[k3p].resize(k_mesh);
    for (int k3 = 0; k3 < k_mesh; k3++) {
      P3 = Pwave_vecwx[0][n3][l3][k3];
      PglbarP[k3p][k3] = glbar_vec[0][l3][k3p][k3] * P3p * P3;
    }
  }

  temp_2 = 0.0;
  for (int k_a = 0; k_a < k_mesh; k_a++) {
    temp_1 = 0.0;
    for (int k_b = 0; k_b < k_mesh; k_b++) {
      temp_1 += PF12[k_b] * PglbarP[k_a][k_b];
    }
    temp_2 += temp_1 * PF12_p[k_a];
  }
  result = temp_2;

  return result;
}

double Wave_Fun::PFKglbar2_cal2(int &n12p, int &n12, int &n3p, int &n3,
                                int &l3p, int &l3, int lambda1, int &lbar) {
  double result = 0.0;
  int lambda2 = 2 - lambda1;

  double temp_1 = 0.0, temp_2 = 0.0;
  vector<double> PF12_p(k_mesh), PF12(k_mesh);
  vector<vector<double>> PglbarP(k_mesh);

  for (int k_a = 0; k_a < k_mesh; k_a++) {
    temp_1 = 0.0;
    temp_2 = 0.0;
    PglbarP[k_a].resize(k_mesh);
    for (int k_b = 0; k_b < k_mesh; k_b++) {
      temp_1 += Pwave_vecwx[0][n12p][0][k_b] * F_reg_vec[k_a][k_b];
      temp_2 += Pwave_vecwx[0][n12][0][k_b] * F_reg_vec[k_a][k_b];

      PglbarP[k_a][k_b] = glbar_vec[1][lbar][k_a][k_b] *
                          Pwave_vecwx[lambda1][n3p][l3p][k_a] *
                          Pwave_vecwx[lambda2][n3][l3][k_b];
    }
    PF12_p[k_a] = temp_1;
    PF12[k_a] = temp_2;
  }

  temp_2 = 0.0;
  for (int k_a = 0; k_a < k_mesh; k_a++) {
    temp_1 = 0.0;
    for (int k_b = 0; k_b < k_mesh; k_b++) {
      temp_1 += PF12[k_b] * PglbarP[k_a][k_b];
    }
    temp_2 += temp_1 * PF12_p[k_a];
  }
  result = temp_2;

  return result;
}

//  Legendre from net .........................................................

void Wave_Fun::legendre_compute_glr(int n, double x[], double w[])

//****************************************************************************80
//
//  Purpose:
//
//    LEGENDRE_COMPUTE_GLR: Legendre quadrature by the Glaser-Liu-Rokhlin
//    method.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    20 October 2009
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Reference:
//
//    Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
//    A fast algorithm for the calculation of the roots of special functions,
//    SIAM Journal on Scientific Computing,
//    Volume 29, Number 4, pages 1420-1438, 2007.
//
//  Parameters:
//
//    Input, int N, the order.
//
//    Output, double X[N], the abscissas.
//
//    Output, double W[N], the weights.
//
{
  int i;
  double p;
  double pp;
  double w_sum;
  //
  //  Get the value and derivative of the N-th Legendre polynomial at 0.
  //
  legendre_compute_glr0(n, &p, &pp);
  //  cout<<"???"<<p<<"\t"<<pp<<endl;
  //
  //  If N is odd, then zero is a root.
  //
  if (n % 2 == 1) {
    x[(n - 1) / 2] = p;
    w[(n - 1) / 2] = pp;
  }
  //
  //  If N is even, we have to call a function to find the first root.
  //
  else {
    legendre_compute_glr2(p, n, &x[n / 2], &w[n / 2]);
  }
  //
  //  Get the complete set of roots and derivatives.
  //
  legendre_compute_glr1(n, x, w);
  //
  //  Compute the W.
  //
  for (i = 0; i < n; i++) {
    w[i] = 2.0 / (1.0 - x[i]) / (1.0 + x[i]) / w[i] / w[i];
  }
  w_sum = 0.0;
  for (i = 0; i < n; i++) {
    w_sum = w_sum + w[i];
  }
  for (i = 0; i < n; i++) {
    w[i] = 2.0 * w[i] / w_sum;
  }
  return;
}

//****************************************************************************80

void Wave_Fun::legendre_compute_glr0(int n, double *p, double *pp)

//****************************************************************************80
//
//  Purpose:
//
//    LEGENDRE_COMPUTE_GLR0 gets a starting value for the fast algorithm.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 October 2009
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Reference:
//
//    Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
//    A fast algorithm for the calculation of the roots of special functions,
//    SIAM Journal on Scientific Computing,
//    Volume 29, Number 4, pages 1420-1438, 2007.
//
//  Parameters:
//
//    Input, int N, the order of the Legendre polynomial.
//
//    Output, double *P, *PP, the value of the N-th Legendre polynomial
//    and its derivative at 0.
//
{
  double dk;
  int k;
  double pm1;
  double pm2;
  double ppm1;
  double ppm2;

  pm2 = 0.0;
  pm1 = 1.0;
  ppm2 = 0.0;
  ppm1 = 0.0;

  for (k = 0; k < n; k++) {
    dk = (double)k;
    *p = -dk * pm2 / (dk + 1.0);
    *pp = ((2.0 * dk + 1.0) * pm1 - dk * ppm2) / (dk + 1.0);
    pm2 = pm1;
    pm1 = *p;
    ppm2 = ppm1;
    ppm1 = *pp;
  }
  return;
}

void Wave_Fun::legendre_compute_glr1(int n, double x[], double w[])

//****************************************************************************80
//
//  Purpose:
//
//    LEGENDRE_COMPUTE_GLR1 gets the complete set of Legendre points and
//    weights.
//
//  Discussion:
//
//    This routine requires that a starting estimate be provided for one
//    root and its derivative.  This information will be stored in entry
//    (N+1)/2 if N is odd, or N/2 if N is even, of X and W.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 October 2009
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Reference:
//
//    Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
//    A fast algorithm for the calculation of the roots of special functions,
//    SIAM Journal on Scientific Computing,
//    Volume 29, Number 4, pages 1420-1438, 2007.
//
//  Parameters:
//
//    Input, int N, the order of the Legendre polynomial.
//
//    Input/output, double X[N].  On input, a starting value
//    has been set in one entry.  On output, the roots of the Legendre
//    polynomial.
//
//    Input/output, double W[N].  On input, a starting value
//    has been set in one entry.  On output, the derivatives of the Legendre
//    polynomial at the zeros.
//
//  Local Parameters:
//
//    Local, int M, the number of terms in the Taylor expansion.
//
{
  double dk;
  double dn;
  double h;
  int j;
  int k;
  int l;
  int m = 30;
  int n2;
  static double pi = 3.141592653589793;
  int s;
  double *u;
  double *up;
  double xp;

  if (n % 2 == 1) {
    n2 = (n - 1) / 2 - 1;
    s = 1;
  } else {
    n2 = n / 2 - 1;
    s = 0;
  }

  u = new double[m + 2];
  up = new double[m + 1];

  dn = (double)n;

  for (j = n2 + 1; j < n - 1; j++) {
    xp = x[j];

    h = rk2_leg(pi / 2.0, -pi / 2.0, xp, n) - xp;

    u[0] = 0.0;
    u[1] = 0.0;
    u[2] = w[j];

    up[0] = 0.0;
    up[1] = u[2];

    for (k = 0; k <= m - 2; k++) {
      dk = (double)k;

      u[k + 3] = (2.0 * xp * (dk + 1.0) * u[k + 2] +
                  (dk * (dk + 1.0) - dn * (dn + 1.0)) * u[k + 1] / (dk + 1.0)) /
                 (1.0 - xp) / (1.0 + xp) / (dk + 2.0);

      up[k + 2] = (dk + 2.0) * u[k + 3];
    }

    for (l = 0; l < 5; l++) {
      h = h - ts_mult(u, h, m) / ts_mult(up, h, m - 1);
    }

    x[j + 1] = xp + h;
    w[j + 1] = ts_mult(up, h, m - 1);
  }

  for (k = 0; k <= n2 + s; k++) {
    x[k] = -x[n - 1 - k];
    w[k] = w[n - 1 - k];
  }
  return;
}

void Wave_Fun::legendre_compute_glr2(double pn0, int n, double *x1, double *d1)

//****************************************************************************80
//
//  Purpose:
//
//    LEGENDRE_COMPUTE_GLR2 finds the first real root.
//
//  Discussion:
//
//    This function is only called if N is even.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 October 2009
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Reference:
//
//    Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
//    A fast algorithm for the calculation of the roots of special functions,
//    SIAM Journal on Scientific Computing,
//    Volume 29, Number 4, pages 1420-1438, 2007.
//
//  Parameters:
//
//    Input, double PN0, the value of the N-th Legendre polynomial
//    at 0.
//
//    Input, int N, the order of the Legendre polynomial.
//
//    Output, double *X1, the first real root.
//
//    Output, double *D1, the derivative at X1.
//
//  Local Parameters:
//
//    Local, int M, the number of terms in the Taylor expansion.
//
{
  double dk;
  double dn;
  int k;
  int l;
  int m = 30;
  static double pi = 3.141592653589793;
  double t;
  double *u;
  double *up;

  t = 0.0;
  *x1 = rk2_leg(t, -pi / 2.0, 0.0, n);

  u = new double[m + 2];
  up = new double[m + 1];

  dn = (double)n;
  //
  //  U[0] and UP[0] are never used.
  //  U[M+1] is set, but not used, and UP[M] is set and not used.
  //  What gives?
  //
  u[0] = 0.0;
  u[1] = pn0;

  up[0] = 0.0;

  for (k = 0; k <= m - 2; k = k + 2) {
    dk = (double)k;

    u[k + 2] = 0.0;
    u[k + 3] = (dk * (dk + 1.0) - dn * (dn + 1.0)) * u[k + 1] / (dk + 1.0) /
               (dk + 2.0);

    up[k + 1] = 0.0;
    up[k + 2] = (dk + 2.0) * u[k + 3];
  }

  for (l = 0; l < 5; l++) {
    *x1 = *x1 - ts_mult(u, *x1, m) / ts_mult(up, *x1, m - 1);
  }
  *d1 = ts_mult(up, *x1, m - 1);

  return;
}

void Wave_Fun::rescale(double a, double b, int n, double x[], double w[])

//****************************************************************************80
//
//  Purpose:
//
//    RESCALE rescales a Legendre quadrature rule from [-1,+1] to [A,B].
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    18 October 2009
//
//  Author:
//
//    Original MATLAB version by Nick Hale.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin,
//    A fast algorithm for the calculation of the roots of special functions,
//    SIAM Journal on Scientific Computing,
//    Volume 29, Number 4, pages 1420-1438, 2007.
//
//  Parameters:
//
//    Input, double A, B, the endpoints of the new interval.
//
//    Input, int N, the order.
//
//    Input/output, double X[N], on input, the abscissas for [-1,+1].
//    On output, the abscissas for [A,B].
//
//    Input/output, double W[N], on input, the weights for [-1,+1].
//    On output, the weights for [A,B].
//
{
  int i;

  for (i = 0; i < n; i++) {
    x[i] = ((a + b) + (b - a) * x[i]) / 2.0;
  }
  for (i = 0; i < n; i++) {
    w[i] = (b - a) * w[i] / 2.0;
  }
  return;
}

//****************************************************************************80

double Wave_Fun::rk2_leg(double t1, double t2, double x, int n)

//****************************************************************************80
//
//  Purpose:
//
//    RK2_LEG advances the value of X(T) using a Runge-Kutta method.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    22 October 2009
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double T1, T2, the range of the integration interval.
//
//    Input, double X, the value of X at T1.
//
//    Input, int N, the number of steps to take.
//
//    Output, double RK2_LEG, the value of X at T2.
//
{
  double f;
  double h;
  int j;
  double k1;
  double k2;
  int m = 10;
  double snn1;
  double t;

  h = (t2 - t1) / (double)m;
  snn1 = sqrt((double)(n * (n + 1)));
  t = t1;

  for (j = 0; j < m; j++) {
    f = (1.0 - x) * (1.0 + x);
    k1 = -h * f / (snn1 * sqrt(f) - 0.5 * x * sin(2.0 * t));
    x = x + k1;

    t = t + h;

    f = (1.0 - x) * (1.0 + x);
    k2 = -h * f / (snn1 * sqrt(f) - 0.5 * x * sin(2.0 * t));
    x = x + 0.5 * (k2 - k1);
  }
  return x;
}

double Wave_Fun::ts_mult(double *u, double h, int n)

//****************************************************************************80
//
//  Purpose:
//
//    TS_MULT evaluates a polynomial.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    17 May 2013
//
//  Author:
//
//    Original C++ version by Nick Hale.
//    This C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, double U[N+1], the polynomial coefficients.
//    U[0] is ignored.
//
//    Input, double H, the polynomial argument.
//
//    Input, int N, the number of terms to compute.
//
//    Output, double TS_MULT, the value of the polynomial.
//
{
  double hk;
  int k;
  double ts;

  ts = 0.0;
  hk = 1.0;
  for (k = 1; k <= n; k++) {
    ts = ts + u[k] * hk;
    hk = hk * h;
  }
  return ts;
}

void Wave_Fun::print() {
  for (int i = 0; i < k_mesh; i++) {
    std::cout << Pwave_x[0][i] << '\t' << Pwave_w[i] << '\t'
              << Pwave_vec[1][0][i] << std::endl;
  }
}
void Wave_Fun::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  for (int i = 0; i < k_mesh; i++) {
    file << Pwave_x[0][i] << '\t' << Pwave_w[i] << '\t' << std::endl;
  }
  file.close();
}
void Wave_Fun::print_Pmesh(std::string filename) {
  std::ofstream file;
  file.open(filename);
  for (int i = 0; i < p_mesh; i++) {
    file << legendreP_x[i] << '\t' << legendreP_w[i] << std::endl;
  }
  file.close();
}
void Wave_Fun::glbarS_init(int lbar, int lbar12, int lbar3,
                           int lmax_legendrePlW) {
  // init PFKglbarS for 2-pi exchange term
  //

  int lmax = std::max(lbar, lbar12);
  lmax = std::max(lmax, lbar3);
  this->lbar_glbarS_max = lmax;

  int lmax_LP = lmax_legendrePlW;
  legendreP_x.resize(p_mesh);
  legendreP_w.resize(p_mesh);
  this->gauleg(-1.0, 1.0, legendreP_x, legendreP_w);

  legendre_PlW.resize(lmax_LP + 2);
  for (int l = 0; l <= lmax_LP + 1; l++) {
    legendre_PlW[l].resize(p_mesh);
    for (int u = 0; u < p_mesh; u++) {
      legendre_PlW[l][u] =
          gsl_sf_legendre_Pl(l, legendreP_x[u]) * legendreP_w[u];
    }
  }

  legendreP_x_L.resize(p_mesh);
  legendreP_w_L.resize(p_mesh);
  long double x_i = -1.0;
  long double x_f = 1.0;
  this->gauleg(x_i, x_f, legendreP_x_L, legendreP_w_L);

  legendre_PlW_L.resize(lmax_LP + 2);
  for (int l = 0; l <= lmax_LP + 1; l++) {
    legendre_PlW_L[l].resize(p_mesh);
    for (int u = 0; u < p_mesh; u++) {
      legendre_PlW_L[l][u] =
          gsl_sf_legendre_Pl(l, legendreP_x_L[u]) * legendreP_w_L[u];
    }
  }
  legendre_PlW_L_B.resize(p_mesh);
  for (int u = 0; u < p_mesh; u++) {
    legendre_PlW_L_B[u].resize(lmax_LP + 2);
    for (int l = 0; l <= lmax_LP + 1; l++) {
      legendre_PlW_L_B[u][l] = legendre_PlW_L[l][u];
    }
  }

  //
  if (WaveCom.rank() == 0) {
    std::ofstream file;
    file.open("data/legendre_check.dat");
    file << "i\t x\t\t x_L\t\t w\t\t w_L" << endl;
    for (int i = 0; i < p_mesh; i++) {
      file << i << "\t " << legendreP_x[i] << "\t " << legendreP_x_L[i];
      file << "\t " << legendreP_w[i] << "\t " << legendreP_w_L[i] << endl;
    }
    file << "\t === PlW === " << endl;
    for (int l = 0; l <= lmax_LP + 1; l++) {
      file << "\t --- l = " << l << " --- " << endl;
      for (int u = 0; u < p_mesh; u++) {
        file << u << "\t " << legendre_PlW[l][u] << "\t "
             << legendre_PlW_L[l][u] << endl;
      }
    }
    file.close();
  }
  // initial lbar_state_g & glbar_vec_g
  // here is smaller that lbar_state & glbar_vec
  int index = 0;
  lbar_struct lbar_state_t;
  lbar_vec_g.resize(lbar + 1);
  lbar_state_g.clear();
  for (int l1 = 0; l1 <= lbar; l1++) {
    lbar_vec_g[l1].resize(lbar12 + 1);
    for (int l2 = 0; l2 <= lbar12; l2++) {
      lbar_vec_g[l1][l2].resize(lbar3 + 1);
      for (int l3 = 0; l3 <= lbar3; l3++) {
        lbar_vec_g[l1][l2][l3] = index;
        lbar_state_t.lbar = l1;
        lbar_state_t.lbar12 = l2;
        lbar_state_t.lbar3 = l3;
        lbar_state_g.push_back(lbar_state_t);
        index++;
      }
    }
  }
  // Slbar_state
  Slbar_vec.resize(4);
  S_lbar_struct Slbar_state_t;
  index = 0;
  for (int S_index = 0; S_index < 4; S_index++) {
    int size_lbar = lbar_state_g.size();
    Slbar_vec[S_index].resize(size_lbar);
    for (int l = 0; l < size_lbar; l++) {
      Slbar_state_t.S_index = S_index;
      Slbar_state_t.lbar_index = l;
      Slbar_state.push_back(Slbar_state_t);
      Slbar_vec[S_index][l] = index;
      index++;
    }
  }
  if (WaveCom.rank() == 0) {
    cout << " ### lbar = " << lbar << "\t lbar12 = " << lbar12
         << "\t lbar3 = " << lbar3
         << "\t lmax_legendrePlW = " << lmax_legendrePlW << endl;
    cout << " ### lbar_state_g.size() = " << lbar_state_g.size() << endl;
  }
  //
  // initial ksq
  int size_kpk = kpk_state.size();
  int size_u = legendreP_x.size();
  if (size_kpk == 0) {
    cout << "Wrong at Wave_Fun::PFKglbarS_init. kpk.size() == 0" << endl;
    exit(0);
  }
  if (size_u == 0) {
    cout << "Wrong at Wave_Fun::PFKglbarS_init. legendreP_x.size() == 0"
         << endl;
    exit(0);
  }
  ksq.resize(size_kpk);
  ksq2.resize(size_kpk);
  for (int i = 0; i < size_kpk; i++) {
    ksq[i].resize(size_u);
    ksq2[i].resize(size_u);
    int n_k = kpk_state[i].n_k;
    int n_kp = kpk_state[i].n_kp;
    double k = Pwave_x[0][n_k];
    double k2 = Pwave_x[1][n_k];
    double kp = Pwave_x[0][n_kp];
    double kp2 = Pwave_x[1][n_kp];
    for (int j = 0; j < size_u; j++) {
      double u = legendreP_x[j];
      double t = k2 + kp2 - 2 * k * kp * u;
      ksq2[i][j] = t;
      ksq[i][j] = sqrt(t);
    }
  }
  ksq_L.resize(size_kpk);
  ksq2_L.resize(size_kpk);
  for (int i = 0; i < size_kpk; i++) {
    ksq_L[i].resize(size_u);
    ksq2_L[i].resize(size_u);
    int n_k = kpk_state[i].n_k;
    int n_kp = kpk_state[i].n_kp;
    double k = Pwave_x[0][n_k];
    double k2 = Pwave_x[1][n_k];
    double kp = Pwave_x[0][n_kp];
    double kp2 = Pwave_x[1][n_kp];
    for (int j = 0; j < size_u; j++) {
      long double u = legendreP_x_L[j];
      long double t = k2 + kp2 - 2 * k * kp * u;
      ksq2_L[i][j] = t;
      ksq_L[i][j] = sqrt(t);
    }
  }
}

double Wave_Fun::glbarS_cal(int s1, int s2, int lbar, int lbar12, int lbar3,
                            int size_u, int kpk_12, int kpk_3) {
  // calculate glbarS
  double result = 0.0;
  // int size_u = legendreP_x.size();
  /*  int lbar=lbar_state[lbar_index].lbar;
    int lbar12=lbar_state[lbar_index].lbar12;
    int lbar3=lbar_state[lbar_index].lbar3;
    int s1,s2;
    switch (S_index) {
      case 0:
        s1 = 0;
        s2 = 0;
        break;
      case 1:
        s1 = 2;
        s2 = 0;
        break;
      case 2:
        s1 = 0;
        s2 = 2;
        break;
      case 3:
        s1 = 2;
        s2 = 2;
        break;
      default:
        cout<<"Wrong at Wave_Fun::glbarS_cal S_index ?"<<endl;
        exit(0);
        break;
    }
  */
  double C1 = 1.0 / 6.0;
  double C2 = 1.0 / std::sqrt(3);
  double C3 = 2.0 / 3.0;

  double t, t1, t2;
  double KP3, KP3_2, KP12, KP12_2, Qp2, Q2;
  double part1, part2, part3, part1_t, part2_t, part3_t;
  // integration over u3 u12 u
  // part3
  part3 = 0.0;
  for (int u3 = 0; u3 < size_u; u3++) {
    // KP3_2=ksq2[kpk_3][u3];
    KP3 = ksq[kpk_3][u3];
    KP3_2 = KP3 * KP3;
    Qp2 = C3 * KP3_2;
    part3_t = 0.0;
    part3_t = 1.0 / (Qp2 + Mass_pi2);
    if (s2 == 0)
      part3_t *= Qp2;
    // part2
    part2 = 0.0;
    t1 = C1 * KP3_2;
    for (int u12 = 0; u12 < size_u; u12++) {
      KP12 = ksq[kpk_12][u12];
      KP12_2 = KP12 * KP12;
      t = KP3 * KP12;
      part2_t = 1.0 / (this->poww(t, lbar));
      // part2_t = 1.0 / (std::pow(t, lbar));
      // part1
      part1 = 0.0;
      t2 = 0.5 * KP12_2;
      for (int u = 0; u < size_u; u++) {
        Q2 = t1 + t2 - C2 * t * legendreP_x[u];
        part1_t = 1.0 / (Q2 + Mass_pi2);
        if (s1 == 0)
          part1_t *= Q2;
        part1 += part1_t * legendre_PlW[lbar][u];
      } // finished the integration of u
      // if (std::abs(part1) < EPS_glbarS)
      //   part1 == 0.0;
      part2 += part2_t * part1 * legendre_PlW[lbar12][u12];
    } // finished the integration of u12
    part3 += part3_t * part2 * legendre_PlW[lbar3][u3];
  } // finished the integration of u3
  // result = part3;
  result = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0 * part3;
  return result;
}

double Wave_Fun::glbarS_cal_for_check(int s1, int s2, int lbar, int lbar12,
                                      int lbar3, int size_u, int kpk_12,
                                      int kpk_3) {
  // calculate glbarS
  double result = 0.0;
  // int size_u = legendreP_x.size();
  /*  int lbar=lbar_state[lbar_index].lbar;
    int lbar12=lbar_state[lbar_index].lbar12;
    int lbar3=lbar_state[lbar_index].lbar3;
    int s1,s2;
    switch (S_index) {
      case 0:
        s1 = 0;
        s2 = 0;
        break;
      case 1:
        s1 = 2;
        s2 = 0;
        break;
      case 2:
        s1 = 0;
        s2 = 2;
        break;
      case 3:
        s1 = 2;
        s2 = 2;
        break;
      default:
        cout<<"Wrong at Wave_Fun::glbarS_cal S_index ?"<<endl;
        exit(0);
        break;
    }
  */
  int p_mesh2 = 50;
  int size_u2 = p_mesh2;
  int lmax_LP = 14;
  vector<double> legendreP_x2, legendreP_w2;
  legendreP_x2.resize(p_mesh2);
  legendreP_w2.resize(p_mesh2);
  this->gauleg(-1.0, 1.0, legendreP_x2, legendreP_w2);

  vector<vector<double>> legendre_PlW2;
  legendre_PlW2.resize(lmax_LP + 2);
  for (int l = 0; l <= lmax_LP + 1; l++) {
    legendre_PlW2[l].resize(p_mesh2);
    for (int u = 0; u < p_mesh2; u++) {
      legendre_PlW2[l][u] =
          gsl_sf_legendre_Pl(l, legendreP_x2[u]) * legendreP_w2[u];
    }
  }

  std::ofstream file;
  file.open("data/glbarS_check.dat");
  file << "s1 : " << s1 << "\t s2 : " << s2 << "\t lbar : " << lbar
       << "\t lbar12 : " << lbar12 << "\t lbar3 : " << lbar3 << endl;
  file << "kpk12 = " << kpk_12 << "\t k12p_n = " << kpk_state[kpk_12].n_kp
       << "\t k12_n = " << kpk_state[kpk_12].n_k << endl;
  file << "kpk3 = " << kpk_3 << "\t k3p_n = " << kpk_state[kpk_3].n_kp
       << "\t k3_n = " << kpk_state[kpk_3].n_k << endl;
  file << "\t ---- k mesh ---- " << endl;
  file << "n\t k\t w" << endl;
  for (int i = 0; i < Pwave_x[0].size(); i++)
    file << i << "\t " << Pwave_x[0][i] << "\t " << Pwave_w[i] << endl;
  file << " \t---- ksq [kpk_3][u] ---- " << endl;
  int kp_n = kpk_state[kpk_3].n_kp;
  file << "i\t u\t sqrt(kp^2 + k^2 - 2*kp*k*u)"
       << "  kp = " << Pwave_x[0][kp_n] << endl;
  for (int i = 0; i < size_u; i++) {
    file << i << "\t " << legendreP_x[i] << "\t " << ksq[kpk_3][i] << endl;
  }
  file << "u3\t u12\t u\t part1_t\t Plu3\t add" << endl;
  double C1 = 1.0 / 6.0;
  double C2 = 1.0 / std::sqrt(3);
  double C3 = 2.0 / 3.0;

  double t, t1, t2;
  double KP3, KP3_2, KP12, KP12_2, Qp2, Q2;
  double part1, part2, part3, part1_t, part2_t, part3_t;
  // integration over u3 u12 u
  // part3
  part3 = 0.0;
  for (int u3 = 0; u3 < size_u; u3++) {
    // KP3_2=ksq2[kpk_3][u3];
    KP3 = ksq[kpk_3][u3];
    KP3_2 = KP3 * KP3;
    Qp2 = C3 * KP3_2;
    part3_t = 0.0;
    part3_t = 1.0 / (Qp2 + Mass_pi2);
    if (s2 == 0)
      part3_t *= Qp2;
    // part2
    part2 = 0.0;
    t1 = C1 * KP3_2;
    for (int u12 = 0; u12 < size_u; u12++) {
      KP12 = ksq[kpk_12][u12];
      KP12_2 = KP12 * KP12;
      t = KP3 * KP12;
      part2_t = 1.0 / (this->poww(t, lbar));
      // part2_t = 1.0 / (std::pow(t, lbar));
      // part1
      part1 = 0.0;
      t2 = 0.5 * KP12_2;
      for (int u = 0; u < size_u2; u++) {
        Q2 = t1 + t2 - C2 * t * legendreP_x2[u];
        part1_t = 1.0 / (Q2 + Mass_pi2);
        if (s1 == 0)
          part1_t *= Q2;
        double add = part1_t * legendre_PlW2[lbar][u];
        part1 += add;
        file << u3 << "\t " << u12 << "\t " << u << "\t\t " << part1_t << "\t "
             << legendre_PlW2[lbar][u] << "\t " << add << endl;
      } // finished the integration of u
      part2 += part2_t * part1 * legendre_PlW[lbar12][u12];
      file << " ------ part2 : part2_t=" << part2_t << "\t * part1=" << part1
           << "\t * Pl12=" << legendre_PlW[lbar12][u12] << "  =   "
           << part2_t * part1 * legendre_PlW[lbar12][u12] << endl;
    } // finished the integration of u12
    part3 += part3_t * part2 * legendre_PlW[lbar3][u3];
    file << " ====== part3 : part3_t=" << part3_t << "\t * part2=" << part2
         << "\t * Pl3=" << legendre_PlW[lbar3][u3] << "  =   "
         << part3_t * part2 * legendre_PlW[lbar3][u3] << endl;
  } // finished the integration of u3
  // result = part3;
  result = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0 * part3;
  file.close();
  cout << " ====== FINISHED CHECK glbarS ======= " << endl;
  // exit(0);

  return result;
}

double Wave_Fun::poww(double &a, int b) {
  double ans = 1.0;
  double base = a;
  while (b != 0) {
    if (b & 1 != 0)
      ans *= base;
    base *= base;
    b >>= 1;
  }
  return ans;
}
void Wave_Fun::glbarS_vec_build() {
  // glbar_vec[S_index][lbar_state][kpk_12][kpk_3]
  // S_index: s1   s2    0.5*s1 + s2
  //         0    0            0
  //         2    0            1
  //         0    2            2
  //         2    2            3
  glbarS_vec.resize(4);
  glbarS2_vec.resize(4);
  // double G = 0.0;
  // int k12, k12p, k3, k3p;
  // int k123_index, k123p_index;
  int s1, s2;
  int size_lbar = lbar_state_g.size();
  int size_kpk = kpk_state.size();
  cout << "Begin to glbarS_vec_build " << endl;
  for (int S_index = 0; S_index < 4; S_index++) {
    switch (S_index) {
    case 0:
      s1 = 0;
      s2 = 0;
      break;
    case 1:
      s1 = 2;
      s2 = 0;
      break;
    case 2:
      s1 = 0;
      s2 = 2;
      break;
    case 3:
      s1 = 2;
      s2 = 2;
      break;
    default:
      cout << "Wrong at Wave_Fun::glbarS_vec_build S_index ?" << endl;
      exit(0);
      break;
    }
    glbarS_vec[S_index].resize(size_lbar);
    glbarS2_vec[S_index].resize(size_lbar);
    for (int l = 0; l < size_lbar; l++) {
      glbarS_vec[S_index][l].resize(size_kpk);
      glbarS2_vec[S_index][l].resize(size_kpk);
      for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
        glbarS_vec[S_index][l][kpk12].resize(size_kpk);
        glbarS2_vec[S_index][l][kpk12].resize(size_kpk);
      }
    }
    //#pragma omp parallel for
    // cout << "S_index : " << S_index << endl;
    Pgs_Bar pgs;
    // pgs.init(size_lbar, 30);
    for (int l = 0; l < size_lbar; l++) {
      // pgs.PrintBar(l);
      // cout << l << endl;
      int lbar = lbar_state_g[l].lbar;
      int lbar12 = lbar_state_g[l].lbar12;
      int lbar3 = lbar_state_g[l].lbar3;
      double C = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0;
#pragma omp parallel for schedule(dynamic, 1)
      for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
#pragma omp parallel for schedule(dynamic, 1)
        for (int kpk3 = 0; kpk3 < size_kpk; kpk3++) {
          // G = C * glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh, kpk12,
          // kpk3);

          int k12p = kpk_state[kpk12].n_kp;
          int k12 = kpk_state[kpk12].n_k;
          int k3p = kpk_state[kpk3].n_kp;
          int k3 = kpk_state[kpk3].n_k;
          int k123p_index = k123_vec[k12p][k3p];
          int k123_index = k123_vec[k12][k3];
          double G =
              glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh, kpk12, kpk3);
          glbarS_vec[S_index][l][kpk12][kpk3] = G;
          // cout<<kpk12<<"\t"<<kpk3<<"\t"<<k123_index<<"\t"<<k123p_index<<endl;
          glbarS2_vec[S_index][l][k123p_index][k123_index] = G;
        }
      }
    }
    // pgs.EndBar();
  }
}

void Wave_Fun::glbarS_vec_build_quick() {
  // glbar_vec[S_index][lbar_state][kpk_12][kpk_3]
  // S_index: s1   s2    0.5*s1 + s2
  //         0    0            0
  //         2    0            1
  //         0    2            2
  //         2    2            3
  const int size_u = p_mesh;
  const int lbar_max_LP = lbar_glbarS_max;
  const int size_lbar = lbar_state_g.size();
  const int size_kpk = kpk_state.size();
  const long double C1 = 1.0 / 6.0;
  const long double C2 = 1.0 / std::sqrt(3);
  const long double C3 = 2.0 / 3.0;
  const long double C_one = 1.0;

  glbarS_vec.resize(4);
  glbarS2_vec.resize(4);
  for (int S_index = 0; S_index < 4; S_index++) {
    glbarS_vec[S_index].resize(size_lbar);
    glbarS2_vec[S_index].resize(size_lbar);
    for (int l = 0; l < size_lbar; l++) {
      glbarS_vec[S_index][l].resize(size_kpk);
      glbarS2_vec[S_index][l].resize(size_kpk);
      for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
        glbarS_vec[S_index][l][kpk12].resize(size_kpk);
        glbarS2_vec[S_index][l][kpk12].resize(size_kpk);
      }
    }
  }

  //
  // Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar_p, Leg_mat_lbar_m;
  // Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar3_p, Leg_mat_lbar3_m;
  // Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar12_p, Leg_mat_lbar12_m;
  Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar(size_lbar, size_u);
  Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar12(size_lbar, size_u);
  Matrix<long double, Dynamic, Dynamic> Leg_mat_lbar3(size_lbar, size_u);

  // Leg_mat_lbar12_p.resize(size_u, size_lbar);
  // Leg_mat_lbar12_m.resize(size_u, size_lbar);
  Leg_mat_lbar12.resize(size_lbar, size_u);
  // Leg_mat_lbar3_p.resize(size_u, size_lbar);
  // Leg_mat_lbar3_m.resize(size_u, size_lbar);
  Leg_mat_lbar3.resize(size_lbar, size_u);
  for (int l = 0; l < size_lbar; l++) {
    int lbar = lbar_state_g[l].lbar;
    int lbar12 = lbar_state_g[l].lbar12;
    int lbar3 = lbar_state_g[l].lbar3;
    for (int u = 0; u < size_u; u++) {
      // Leg_mat_lbar_p(l, u) = 0.0;
      // Leg_mat_lbar_m(l, u) = 0.0;
      long double L_lbar = legendre_PlW_L[lbar][u];
      // Leg_mat_lbar(l, u) = L_lbar;
      // if (L_lbar < 0) {
      //   Leg_mat_lbar_m(l, u) = L_lbar;
      // } else {
      //   Leg_mat_lbar_p(l, u) = L_lbar;
      // }
      Leg_mat_lbar(l, u) = legendre_PlW_L[lbar][u];
      Leg_mat_lbar12(l, u) = legendre_PlW_L[lbar12][u];
      Leg_mat_lbar3(l, u) = legendre_PlW_L[lbar3][u];
    }
  }
  //
  Matrix<long double, Dynamic, 1> One_mat_u(size_u);
  for (int i = 0; i < size_u; i++)
    One_mat_u(i) = 1.0;
  //
  Matrix<long double, Dynamic, 1> G_mat(size_lbar);
  for (int l = 0; l < size_lbar; l++) {
    const int lbar = lbar_state_g[l].lbar;
    const int lbar12 = lbar_state_g[l].lbar12;
    const int lbar3 = lbar_state_g[l].lbar3;
    G_mat(l) = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0;
  }
  //

  cout << "Begin to glbarS_cal_Init " << endl;
  // this->glbarS_cal_Init();
  cout << "Begin to glbarS_vec_build " << endl;
  Pgs_Bar pgs;
  // pgs.init(size_lbar, 30);
#pragma omp parallel for schedule(dynamic, 1)
  for (int kpk3 = 0; kpk3 < size_kpk; kpk3++) {
    vector<long double> &ksq_3 = ksq_L[kpk3];
#pragma omp parallel for schedule(dynamic, 1)
    for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
      vector<long double> &ksq_12 = ksq_L[kpk12];
      const int k12p = kpk_state[kpk12].n_kp;
      const int k12 = kpk_state[kpk12].n_k;
      const int k3p = kpk_state[kpk3].n_kp;
      const int k3 = kpk_state[kpk3].n_k;
      const int k123p_index = k123_vec[k12p][k3p];
      const int k123_index = k123_vec[k12][k3];

      Matrix<long double, Dynamic, 1> part3t_mat_s2_0(size_u);
      Matrix<long double, Dynamic, 1> part3t_mat_s2_2(size_u);
      Matrix<long double, Dynamic, Dynamic> part3t_mat_s1_0(size_lbar, size_u);
      Matrix<long double, Dynamic, Dynamic> part3t_mat_s1_2(size_lbar, size_u);
      Matrix<long double, Dynamic, 1> part3_mat_s1_0_s2_0(size_lbar);
      Matrix<long double, Dynamic, 1> part3_mat_s1_0_s2_2(size_lbar);
      Matrix<long double, Dynamic, 1> part3_mat_s1_2_s2_0(size_lbar);
      Matrix<long double, Dynamic, 1> part3_mat_s1_2_s2_2(size_lbar);

      //#pragma omp parallel for reduction(+ : part3)
      for (int u3 = 0; u3 < size_u; u3++) {

        const long double KP3 = ksq_3[u3];
        const long double KP3_2 = KP3 * KP3;
        const long double A = C1 * KP3_2;

        Matrix<long double, Dynamic, 1> part2_mat_s1_0(size_lbar);
        // Matrix<long double, Dynamic, 1> part2_s1_0_m(size_lbar);
        Matrix<long double, Dynamic, 1> part2_mat_s1_2(size_lbar);
        // Matrix<long double, Dynamic, 1> part2_s1_2_m(size_lbar);
        Matrix<long double, Dynamic, Dynamic> part2t_mat_s1_0(size_lbar,
                                                              size_u);
        Matrix<long double, Dynamic, Dynamic> part2t_mat_s1_2(size_lbar,
                                                              size_u);
        // vector<vector<long double>> g_P1_3 = g_P1[u3];
        for (int u12 = 0; u12 < size_u; u12++) {
          // part2_t = glbarS_p2_u3u12[lbar][kpk_3][kpk_12][u3][u12];
          // part2_t = g_P2[u3][u12];
          const long double KP12 = ksq_12[u12];
          const long double KP12_2 = KP12 * KP12;
          const long double B = 0.5 * KP12_2;
          const long double t_12_3 = KP12 * KP3;
          const long double C = C2 * t_12_3;

          Matrix<long double, Dynamic, 1> part1t_mat_s1_0(size_u);
          Matrix<long double, Dynamic, 1> part1t_mat_s1_2(size_u);
          Matrix<long double, Dynamic, 1> part1_mat_s1_0(size_lbar);
          Matrix<long double, Dynamic, 1> part1_mat_s1_2(size_lbar);
          for (int u = 0; u < size_u; u++) {
            const long double Q2 = A + B - C * legendreP_x_L[u];
            const long double dom = Q2 + Mass_pi2;
            const long double part1_t = 1.0 / dom;
            // if (s1 == 0)
            //   part1_t *= Q2;
            part1t_mat_s1_0(u) = part1_t * Q2;
            part1t_mat_s1_2(u) = part1_t;
          }
          part1_mat_s1_0 = Leg_mat_lbar * part1t_mat_s1_0; // +
          // Part1_mat_s1_0 = Part1_mat_s1_0 + Leg_mat_lbar_m * Part1t_mat_s1_0;
          part1_mat_s1_2 = Leg_mat_lbar * part1t_mat_s1_2; // +
          // Part1_mat_s1_2 = Part1_mat_s1_2 + Leg_mat_lbar_m * Part1t_mat_s1_2;
          Matrix<long double, Dynamic, 1> g_P2(lbar_max_LP + 1);
          g_P2(0) = 1.0;
          const long double T = 1.0 / t_12_3;
          for (int i = 1; i <= lbar_max_LP; i++) {
            g_P2(i) = g_P2(i - 1) * T;
          }
          for (int l = 0; l < size_lbar; l++) {
            const int lbar = lbar_state_g[l].lbar;
            const long double gPL = g_P2(lbar) * Leg_mat_lbar12(l, u12);
            part2t_mat_s1_0(l, u12) = gPL * part1_mat_s1_0(l);
            part2t_mat_s1_2(l, u12) = gPL * part1_mat_s1_2(l);
          }
        } // finished the integration of u12

        part2_mat_s1_0 = part2t_mat_s1_0 * One_mat_u;
        part2_mat_s1_2 = part2t_mat_s1_2 * One_mat_u;
        // part2_s1_0 = part2_s1_0 + part2_s1_0_m;
        // part2_s1_2 = part2_s1_2 + part2_s1_2_m;
        // part2 = part2_m + part2_p;
        // long double part3_t = g_P3[u3];
        const long double Q3 = C3 * KP3_2;
        const long double part3_t_s2_2 = 1.0 / (Q3 + Mass_pi2);
        part3t_mat_s2_2(u3) = part3_t_s2_2;
        const long double part3_t_s2_0 = part3_t_s2_2 * Q3;
        part3t_mat_s2_0(u3) = part3_t_s2_0;
        // vector<long double> &Leg_P3 = legendre_PlW_L_B[u3];
        for (int l = 0; l < size_lbar; l++) {
          part3t_mat_s1_0(l, u3) = part2_mat_s1_0(l) * Leg_mat_lbar3(l, u3);
          part3t_mat_s1_2(l, u3) = part2_mat_s1_2(l) * Leg_mat_lbar3(l, u3);
        }
      } // finished the integration of u3
      part3_mat_s1_0_s2_0 = part3t_mat_s1_0 * part3t_mat_s2_0;
      part3_mat_s1_2_s2_0 = part3t_mat_s1_2 * part3t_mat_s2_0;
      part3_mat_s1_0_s2_2 = part3t_mat_s1_0 * part3t_mat_s2_2;
      part3_mat_s1_2_s2_2 = part3t_mat_s1_2 * part3t_mat_s2_2;

      for (int l = 0; l < size_lbar; l++) {
        // const int lbar = lbar_state_g[l].lbar;
        // const int lbar12 = lbar_state_g[l].lbar12;
        // const int lbar3 = lbar_state_g[l].lbar3;
        // double G_c = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1)
        // / 8.0;
        const long double G_c = G_mat(l);
        const double G_00 = G_c * part3_mat_s1_0_s2_0(l);
        const double G_02 = G_c * part3_mat_s1_0_s2_2(l);
        const double G_20 = G_c * part3_mat_s1_2_s2_0(l);
        const double G_22 = G_c * part3_mat_s1_2_s2_2(l);
        glbarS_vec[0][l][kpk12][kpk3] = G_00;
        glbarS_vec[1][l][kpk12][kpk3] = G_20;
        glbarS_vec[2][l][kpk12][kpk3] = G_02;
        glbarS_vec[3][l][kpk12][kpk3] = G_22;
        // cout<<kpk12<<"\t"<<kpk3<<"\t"<<k123_index<<"\t"<<k123p_index<<endl;
        glbarS2_vec[0][l][k123p_index][k123_index] = G_00;
        glbarS2_vec[1][l][k123p_index][k123_index] = G_20;
        glbarS2_vec[2][l][k123p_index][k123_index] = G_02;
        glbarS2_vec[3][l][k123p_index][k123_index] = G_22;
      }
    }
  }

  // pgs.EndBar();
}

void Wave_Fun::glbarS_vec_cal(int &S_index, int &lbar, int &lbar12, int &lbar3,
                              vector<vector<double>> &glbarS2_v) {
  // glbar_vec[S_index][kpk_12][kpk_3]
  // S_index: s1   s2    0.5*s1 + s2
  //         0    0            0
  //         2    0            1
  //         0    2            2
  //         2    2            3
  // glbarS2_v.resize(4);
  // double G = 0.0;
  // int k12, k12p, k3, k3p;
  // int k123_index, k123p_index;

  int size_kpk = kpk_state.size();
  // cout << "Begin to glbarS_vec_cal " << endl;
  // for (int S_index = 0; S_index < 4; S_index++) {
  int s1 = s1_vec[S_index];
  int s2 = s2_vec[S_index];
  glbarS2_v.clear();
  glbarS2_v.resize(size_kpk);
  for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
    glbarS2_v[kpk12].resize(size_kpk);
  }
// cout << l << endl;
// int lbar = lbar_state[l].lbar;
// int lbar12 = lbar_state[l].lbar12;
// int lbar3 = lbar_state[l].lbar3;
// double C = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0;
#pragma omp parallel for schedule(dynamic, 1)
  for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
#pragma omp parallel for schedule(dynamic, 1)
    for (int kpk3 = 0; kpk3 < size_kpk; kpk3++) {
      // G = C * glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh, kpk12,
      // kpk3);

      int k12p = kpk_state[kpk12].n_kp;
      int k12 = kpk_state[kpk12].n_k;
      int k3p = kpk_state[kpk3].n_kp;
      int k3 = kpk_state[kpk3].n_k;
      int k123p_index = k123_vec[k12p][k3p];
      int k123_index = k123_vec[k12][k3];
      double G = glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh, kpk12, kpk3);
      glbarS2_v[k123p_index][k123_index] = G;
    }
  }
}

void Wave_Fun::glbarS_vec_build_MPI_II() {
  // glbar_vec[S_index][lbar_state][kpk_12][kpk_3]
  // S_index: s1   s2    0.5*s1 + s2
  //         0    0            0
  //         2    0            1
  //         0    2            2
  //         2    2            3

  // glbarS_vec.resize(4);
  // glbarS2_vec.resize(4);
  // double G = 0.0;
  // int k12, k12p, k3, k3p;
  // int k123_index, k123p_index;
  int s1, s2;
  int size_lbar = lbar_state_g.size();
  int size_kpk = kpk_state.size();

  mpi::environment env_wave(mpi::threading::multiple);
  int rank = WaveCom.rank();
  int numprocs = WaveCom.size();
  std::string processor_name = "test"; // env_main.processor_name();

  int sl_size = Slbar_state.size();
  this->mpi_size = WaveCom.size();

  // cout << "\t\t sl_size in glbarS = " << sl_size << endl;

  if (WaveCom.rank() == 0) {
    cout << " sl_size = " << sl_size << "\t mpi_size = " << mpi_size << endl;
  }

  int mpi_layer_max = std::ceil(sl_size * 1.0 / mpi_size);

  vector<int> mpi_layer_vec(mpi_size); // for scater

  for (int mpi_num = 0; mpi_num < mpi_size; mpi_num++) {
    int sl = mpi_num + (mpi_layer_max - 1) * mpi_size;
    if (sl < sl_size) {
      mpi_layer_vec[mpi_num] = mpi_layer_max;
    } else {
      mpi_layer_vec[mpi_num] = mpi_layer_max - 1;
    }
  }

  // int mpi_layer; Don't
  scatter(WaveCom, mpi_layer_vec, mpi_layer, 0);
  // mpi_layer = mpi_layer_vec[rank];
  // vector<vector<vector<double>>> glbarS_mpi(mpi_layer);
  // vector<vector<vector<double>>> glbarS2_mpi(size_kpk);
  glbarS2_mpi.resize(mpi_layer);

  for (int lay = 0; lay < mpi_layer; lay++) {
    // for (int sl = 0; sl < sl_size; sl_size++) {
    glbarS2_mpi[lay].resize(size_kpk);
    for (int i = 0; i < size_kpk; i++) {
      // glbarS_mpi[lay][i].resize(size_kpk);
      glbarS2_mpi[lay][i].resize(size_kpk);
    }

    int sl = lay * mpi_size + WaveCom.rank();
    int S_index = Slbar_state[sl].S_index;
    int l = Slbar_state[sl].lbar_index;
    int s1 = s1_vec[S_index];
    int s2 = s2_vec[S_index];

    int lbar = lbar_state_g[l].lbar;
    int lbar12 = lbar_state_g[l].lbar12;
    int lbar3 = lbar_state_g[l].lbar3;
// double C = (2 * lbar + 1) * (2 * lbar3 + 1) * (2 * lbar12 + 1) / 8.0;
// omp_set_num_threads(8);
#pragma omp parallel for schedule(dynamic, 1) // num_threads(10)
    for (int kpk12 = 0; kpk12 < size_kpk; kpk12++) {
#pragma omp parallel for schedule(dynamic, 1) // num_threads(10)
      for (int kpk3 = 0; kpk3 < size_kpk; kpk3++) {
        // G = C * glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh,
        // kpk12, kpk3);

        // glbarS_vec[S_index][l][kpk12][kpk3] = G;
        int k12p = kpk_state[kpk12].n_kp;
        int k12 = kpk_state[kpk12].n_k;
        int k3p = kpk_state[kpk3].n_kp;
        int k3 = kpk_state[kpk3].n_k;
        int k123p_index = k123_vec[k12p][k3p];
        int k123_index = k123_vec[k12][k3];
        // glbarS2_vec[S_index][l][k123p_index][k123_index] = G;
        double G = glbarS_cal(s1, s2, lbar, lbar12, lbar3, p_mesh, kpk12, kpk3);
        // glbarS_mpi[lay][kpk12][kpk3] = G;
        glbarS2_mpi[lay][k123p_index][k123_index] = G;
      }
    }
  }
  if (WaveCom.rank() == 0) {
    vector<vector<vector<vector<double>>>> glbarS2_Gather;
    gather(WaveCom, glbarS2_mpi, glbarS2_Gather, 0);
    cout << "~~~glbarS Gather finished ~~~" << endl;
    glbarS2_vec.resize(4);
    for (int s = 0; s < 4; s++) {
      glbarS2_vec[s].resize(lbar_state_g.size());
    }
    for (int n = 0; n < glbarS2_Gather.size(); n++) {
      cout << " n = " << n
           << "\t glbarS2_Gather size = " << glbarS2_Gather[n].size() << endl;
    }
    for (int lay = 0; lay < mpi_layer_max; lay++) {
      for (int mpi_num = 0; mpi_num < mpi_size; mpi_num++) {
        int sl = lay * mpi_size + mpi_num;
        if (sl < sl_size) {
          int S_index = Slbar_state[sl].S_index;
          int l = Slbar_state[sl].lbar_index;
          // cout << "S_index : " << S_index << "\tl = " << l << endl;
          //  glbarS2_vec[S_index].push_back(glbarS2_Gather[mpi_num][lay]);
          glbarS2_vec[S_index][l] = glbarS2_Gather[mpi_num][lay];
        }
      }
    }
    // cout << "\t TEST : rank = 1 lay = 0 : " << glbarS2_Gather[1][0][0][0]
    //     << endl;
  } else {
    cout << "\t rank = " << WaveCom.rank() << "\t glbarS Finished";
    cout << "\t on \t" << env_wave.processor_name() << endl;
    gather(WaveCom, glbarS2_mpi, 0);
    if (WaveCom.rank() == 1) {
      // double G = glbarS_cal(0, 0, 0, 0, 1, p_mesh, 0, 0);
      // cout << "\t TEST : rank = 1 lay = 0 : " << glbarS2_mpi[0][0][0]
      //     << "\t G = " << G << endl;
    }
  }
  vector<int> mpi_run_flag_gather(WaveCom.size());
  int mpi_run_flag = 1;
  all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);
  // Start to spread glbarS2_vec
  /*
    int trans_num = 0;
    while (std::pow(2, trans_num) < numprocs) {
      trans_num++;
    }
    trans_num--;

    vector<vector<int>> send_node_vec(trans_num + 1);
    vector<vector<int>> recv_node_vec(trans_num + 1);
    send_node_vec[0].resize(1);
    recv_node_vec[0].resize(1);
    send_node_vec[0][0] = 0;
    recv_node_vec[0][0] = numprocs / 2;
    if (rank == 0) {
      cout << "trans_num = " << trans_num << endl;
    }
    for (int i = 1; i <= trans_num; i++) {
      int size_ii = send_node_vec[i - 1].size();
      //  cout << "size_ii = " << size_ii << endl;
      for (int j = 0; j < size_ii; j++) {
        int a = send_node_vec[i - 1][j];
        int b = recv_node_vec[i - 1][j];
        if (a < 0 || b < 0) {
          cout << "Wrong happened at glbarS_vec_build_MPI_III " << endl;
          exit(0);
        }
        send_node_vec[i].push_back(a);
        send_node_vec[i].push_back(b);
      }
      int size_j = send_node_vec[i].size();
      //  cout << "size_j = " << size_j << endl;
      recv_node_vec[i].resize(size_j);
      int index_last = size_j - 1;
      for (int j = 0; j < index_last; j++) {
        int diff = send_node_vec[i][j + 1] - send_node_vec[i][j];
        if (diff < 2) {
          recv_node_vec[i][j] = -1;
        } else {
          recv_node_vec[i][j] = send_node_vec[i][j] + diff / 2;
        }
      }

      int diff_last = numprocs - 1 - send_node_vec[i][index_last];
      if (send_node_vec[i][index_last] > numprocs - 1) {
        recv_node_vec[i][index_last] = -1;
      } else {
        recv_node_vec[i][index_last] =
            send_node_vec[i][index_last] + (diff_last + 1) / 2;
      } /* - */

  /*for (int j = 0; j < send_node_vec[i].size(); j++) {
    cout << send_node_vec[i][j] << "\t";
  }
  cout << "\t\t\t ";
  for (int j = 0; j < recv_node_vec[i].size(); j++) {
    cout << recv_node_vec[i][j] << "\t";
  }
  cout << endl;/* - *//*
  }
  /* - *//*
  all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);

  if (rank == 0) {
    for (int i = 0; i <= trans_num; i++) {
      for (int j = 0; j < send_node_vec[i].size(); j++) {
        cout << send_node_vec[i][j] << "\t";
      }
      cout << "\t\t\t ";
      for (int j = 0; j < recv_node_vec[i].size(); j++) {
        cout << recv_node_vec[i][j] << "\t";
      }
      cout << endl;
    }
  } /* - *//*
  int rank_send_0 = send_node_vec[0][0];
  int rank_recv_0 = recv_node_vec[0][0];
  if (rank == rank_send_0) {
    WaveCom.send(rank_recv_0, 0, glbarS2_vec);
    cout << "\t " << rank_send_0 << " ->  " << rank_recv_0;
  }
  if (rank == rank_recv_0) {
    WaveCom.recv(rank_send_0, 0, glbarS2_vec);
    cout << endl;
  }
  all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);
  // for (int i = 1; i <= trans_num; i++) {
  for (int i = 1; i <= 2; i++) {
    // cout<<" i = "<<endl;
    for (int j = 0; j < send_node_vec[i].size(); j++) {
      int rank_send = send_node_vec[i][j];
      int rank_recv = recv_node_vec[i][j];
      if (rank_send < 0 || rank_recv < 0)
        continue;
      if (rank == rank_send) {
        WaveCom.isend(rank_recv, i, glbarS2_vec);
        cout << "\t " << rank_send << " ->  " << rank_recv;
      }
      if (rank == rank_recv) {
        WaveCom.irecv(rank_send, i, glbarS2_vec);
        cout << endl;
      }
      all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);
    }
    all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);
  }
  all_gather(WaveCom, mpi_run_flag, mpi_run_flag_gather);/* - */
  vector<vector<vector<vector<double>>>> glbarS2_vec_Part1, glbarS2_vec_Part2,
      glbarS2_vec_Part3, glbarS2_vec_Part4;
  if (WaveCom.rank() == 0) {
    int size = glbarS2_vec.size();
    // int size_half = size / 2;
    int size_1_4 = size / 4;
    int size_2_4 = size / 2;
    int size_3_4 = size_1_4 + size_2_4;
    for (int i = 0; i < size_1_4; i++) {
      glbarS2_vec_Part1.push_back(glbarS2_vec[i]);
    }
    for (int i = size_1_4; i < size_2_4; i++) {
      glbarS2_vec_Part2.push_back(glbarS2_vec[i]);
    }
    for (int i = size_2_4; i < size_3_4; i++) {
      glbarS2_vec_Part3.push_back(glbarS2_vec[i]);
    }
    for (int i = size_3_4; i < size; i++) {
      glbarS2_vec_Part4.push_back(glbarS2_vec[i]);
    }
    cout << "Rank 0 glbarS2_vec.size() : " << size << endl;
  }
  // boost::mpi::broadcast(WaveCom, glbarS2_vec, 0);
  // if (WaveCom.rank() == 0)
  cout << "Start broadcast(WaveCom, glbarS2_vec_Part1, 0);" << endl;
  boost::mpi::broadcast(WaveCom, glbarS2_vec_Part1, 0);
  cout << "Start broadcast(WaveCom, glbarS2_vec_Part2, 0);" << endl;
  boost::mpi::broadcast(WaveCom, glbarS2_vec_Part2, 0);
  cout << "Start broadcast(WaveCom, glbarS2_vec_Part3, 0);" << endl;
  boost::mpi::broadcast(WaveCom, glbarS2_vec_Part3, 0);
  cout << "Start broadcast(WaveCom, glbarS2_vec_Part4, 0);" << endl;
  boost::mpi::broadcast(WaveCom, glbarS2_vec_Part4, 0);
  if (WaveCom.rank() != 0) {
    for (int i = 0; i < glbarS2_vec_Part1.size(); i++) {
      glbarS2_vec.push_back(glbarS2_vec_Part1[i]);
    }
    for (int i = 0; i < glbarS2_vec_Part2.size(); i++) {
      glbarS2_vec.push_back(glbarS2_vec_Part2[i]);
    }
    for (int i = 0; i < glbarS2_vec_Part3.size(); i++) {
      glbarS2_vec.push_back(glbarS2_vec_Part3[i]);
    }
    for (int i = 0; i < glbarS2_vec_Part4.size(); i++) {
      glbarS2_vec.push_back(glbarS2_vec_Part4[i]);
    }
  }
  glbarS2_vec_Part1.clear();
  glbarS2_vec_Part2.clear();
  glbarS2_vec_Part3.clear();
  glbarS2_vec_Part4.clear();
  cout << "glbarS2_vec size = " << glbarS2_vec.size();
  cout << "\t on rank = " << WaveCom.rank() << endl;
  // exit(0);
}

void Wave_Fun::PFKglbarS_init(int n_max, int lbar_max, int Emax_F, int Emax_I,
                              vector<vector<int>> E_JC,
                              vector<l_struct_wave> l_state_loc,
                              vector<vector<vector<vector<int>>>> l_vec_loc,
                              vector<lbar_struct> lbar_state_loc,
                              vector<vector<vector<int>>> lbar_vec_loc,
                              vector<vector<int>> Llbar_Ic_Reduced_loc,
                              vector<vector<int>> Llbar_Ic_Red_vec_loc) {

  int check_node = -1;
  // E_max used as a constrain to nlpnl_state
  this->lbar_PFglbarS_max = lbar_max;
  this->l_state = l_state_loc;
  this->l_vec = l_vec_loc;
  this->lbar_state = lbar_state_loc;
  this->lbar_vec = lbar_vec_loc;
  this->Llbar_Ic_Reduced = Llbar_Ic_Reduced_loc;
  this->Llbar_Ic_Red_vec = Llbar_Ic_Red_vec_loc;
  // ---------- initial nl_vec and nl_state for PFKglbar_vec ---------- //
  cout << "  initial nl_vec and nl_state for PFKglbar_vec  " << endl;
  int l_max = 0;
  int l_min = 1000;
  int size_l = l_state_loc.size();
  for (int i = 0; i < size_l; i++) {
    int l12p = l_state_loc[i].l12p;
    int l3p = l_state_loc[i].l3p;
    int l12 = l_state_loc[i].l12;
    int l3 = l_state_loc[i].l3;
    if (l12p > l_max)
      l_max = l12p;
    if (l3p > l_max)
      l_max = l3p;
    if (l12 > l_max)
      l_max = l12;
    if (l3 > l_max)
      l_max = l3;
    //
    if (l12p < l_min)
      l_min = l12p;
    if (l3p < l_min)
      l_min = l3p;
    if (l12 < l_min)
      l_min = l12;
    if (l3 < l_min)
      l_min = l3;
  }
  //
  if (WaveCom.rank() == check_node) {
    cout << "\t \t @@ l_min = " << l_min << "\t l_max = " << l_max
         << "\t n_max = " << n_max << endl;
  }
  nl_state.clear();
  nl_vec.clear();

  int index = 0;
  nl_vec.resize(n_max + 1);
  for (int n12 = 0; n12 <= n_max; n12++) {
    nl_vec[n12].resize(n_max + 1);
    for (int n3 = 0; n3 <= n_max; n3++) {
      nl_vec[n12][n3].resize(l_max + 1);
      for (int l12 = 0; l12 <= l_max; l12++) {
        nl_vec[n12][n3][l12].resize(l_max + 1);
        for (int l3 = 0; l3 <= l_max; l3++) {
          nl_vec[n12][n3][l12][l3] = -1;
          if (l12 < l_min || l12 > l_max)
            continue;
          if (l3 < l_min || l3 > l_max)
            continue;
          nl_struct nl_state_t;
          nl_state_t.n12 = n12;
          nl_state_t.n3 = n3;
          nl_state_t.l12 = l12;
          nl_state_t.l3 = l3;
          // nl_vec[n12][n3][l12][l3] = index;
          nl_state.push_back(nl_state_t);
          index++;
        }
      }
    }
  }

  //

  int size_nl = nl_state.size();
  vector<int> nl2l_flag_vec(size_nl, -1);
  for (int i = 0; i < size_nl; i++) {
    // for (int i = 0; i < 2; i++) {
    int l12 = nl_state[i].l12;
    int l3 = nl_state[i].l3;
    int flag = -1;
    for (int j = 0; j < l_state_loc.size(); j++) {
      int l12_f = l_state_loc[j].l12p;
      int l3_f = l_state_loc[j].l3p;
      int l12_i = l_state_loc[j].l12;
      int l3_i = l_state_loc[j].l3;

      int flag_f = -1;
      if (l12 == l12_f && l3 == l3_f) {
        flag_f = 1;
      }
      int flag_i = -1;
      if (l12 == l12_i && l3 == l3_i)
        flag_i = 1;

      if (flag_f < 0 && flag_i < 0)
        continue;

      if (flag == -1) {
        if (flag_f > 0 && flag_i < 0) {
          flag = 1;
        }
        if (flag_f < 0 && flag_i > 0) {
          flag = 2;
        }
        if (flag_f > 0 && flag_i > 0) {
          flag = 3;
        }
      }
      if (flag == 1) {
        if (flag_i > 0) {
          flag = 3;
        }
      }
      if (flag == 2) {
        if (flag_f > 0) {
          flag = 3;
        }
      }
      if (flag == 3)
        break;
      // if (flag > 0) {
      //   cout << "\t \t j = " << j << "\t\t  flag_f = " << flag_f
      //        << "\t flag_i = " << flag_i << endl;
      //   break;
      // }
    }
    // if (flag < 0) {
    //   cout << "Wrong happened at Wave_Fun::PFKglbarS_init I on rank : "
    //        << WaveCom.rank() << endl;
    //   exit(0);
    // }
    nl2l_flag_vec[i] = flag;
  }

  if (WaveCom.rank() == check_node) {
    cout << "i\t l12p\t l3p\t l12\t l3" << endl;
    for (int i = 0; i < l_state_loc.size(); i++) {
      cout << i << "\t ";
      cout << l_state_loc[i].l12p << "\t ";
      cout << l_state_loc[i].l3p << "\t ";
      cout << l_state_loc[i].l12 << "\t ";
      cout << l_state_loc[i].l3 << endl;
    }

    cout << "i\t n12\t n3\t l12\t l3\t \t\t nl->l flag" << endl;
    for (int i = 0; i < nl_state.size(); i++) {
      cout << i << "\t " << nl_state[i].n12 << "\t " << nl_state[i].n3 << "\t "
           << nl_state[i].l12 << "\t " << nl_state[i].l3 << "\t"
           << nl2l_flag_vec[i] << endl;
    }
  }
  // exit(0);
  vector<nl_struct> nl_s;
  vector<int> nl2l_flag_t;
  nl_s = nl_state;
  nl2l_flag_t = nl2l_flag_vec;
  nl_state.clear();
  nl2l_flag_vec.clear();
  index = 0;
  for (int i = 0; i < size_nl; i++) {
    if (nl2l_flag_t[i] > 0) {
      nl_struct nl_state_t = nl_s[i];
      int n12 = nl_state_t.n12;
      int n3 = nl_state_t.n3;
      int l12 = nl_state_t.l12;
      int l3 = nl_state_t.l3;
      nl_vec[n12][n3][l12][l3] = index;
      nl_state.push_back(nl_state_t);
      nl2l_flag_vec.push_back(nl2l_flag_t[i]);
      index++;
    }
  }
  nl_s.clear();
  if (WaveCom.rank() == check_node) {
    cout << "i\t l12\t l3\t \t\t nl->l flag" << endl;
    for (int i = 0; i < nl_state.size(); i++) {
      cout << i << "\t " << nl_state[i].l12 << "\t " << nl_state[i].l3 << "\t "
           << nl2l_flag_vec[i] << endl;
    }
  }
  // ---------------------------------------------------------------------//

  // -------------------------- calculate nlpnl_state --------------------- //
  //cout << "   calculate nlpnl_state  " << endl;
  nlpnl_struct nlpnl_state_t;
  int index_nlpnl = 0;
  int nl_size = nl_state.size();
  nlpnl_2_lpl.clear();
  nlpnl_state.clear();
  nlpnl_vec.resize(nl_size);
  for (int i = 0; i < nl_size; i++) {
    nlpnl_vec[i].resize(nl_size);
    for (int j = 0; j < nl_size; j++) {
      nlpnl_vec[i][j] = -1;
    }
    int flag_i = -1;
    if (nl2l_flag_vec[i] == 1 || nl2l_flag_vec[i] == 3)
      flag_i = 1;
    if (flag_i < 0)
      continue;

    int l12p = nl_state[i].l12;
    int l3p = nl_state[i].l3;
    int Par_p = this->phase(l12p + l3p);
    int E_p = 2 * (nl_state[i].n12 + nl_state[i].n3) + l12p + l3p;
    vector<int> E_p_pos;
    for (int k1 = 0; k1 < E_JC.size(); k1++) {
      if (E_p == E_JC[k1][0]) {
        E_p_pos.push_back(k1);
      }
    }
    if (E_p_pos.size() == 0)
      continue;
    for (int j = 0; j < nl_size; j++) {
      int flag_j = -1;
      if (nl2l_flag_vec[j] == 2 || nl2l_flag_vec[j] == 3)
        flag_j = 1;
      if (flag_j < 0)
        continue;
      //
      int l12 = nl_state[j].l12;
      int l3 = nl_state[j].l3;
      int index_l = l_vec_loc[l12p][l12][l3p][l3];
      if (index_l < 0)
        continue;
      int Par = this->phase(nl_state[j].l12 + nl_state[j].l3);
      //
      int E = 2 * (nl_state[j].n12 + nl_state[j].n3) + nl_state[j].l12 +
              nl_state[j].l3;
      int E_JC_flag = -1;
      for (int k2 = 0; k2 < E_p_pos.size(); k2++) {
        int a = E_p_pos[k2];
        if (E == E_JC[a][1]) {
          E_JC_flag = 1;
          break;
        }
      }

      // here should not have E_p == E
      if (Par_p == Par) {
        if (E_p <= Emax_F && E <= Emax_I && E_JC_flag == 1) {
          nlpnl_state_t.nl_p = i;
          nlpnl_state_t.nl = j;
          nlpnl_state.push_back(nlpnl_state_t);
          nlpnl_vec[i][j] = index_nlpnl;
          nlpnl_2_lpl.push_back(index_l);
          index_nlpnl++;
        }
      }
    }
  } //
  // cout << "Rank : " << WaveCom.rank() << "\t nlpnl_state.size() = "
  //      << nlpnl_state.size() << endl;
  if (WaveCom.rank() == check_node) {
    int size_t = nlpnl_state.size();
    for (int i = 0; i < size_t; i++) {
      cout << i << "\t " << endl;
      int a = nlpnl_state[i].nl_p;
      int b = nlpnl_state[i].nl;
      int E_p = 2 * (nl_state[a].n12 + nl_state[a].n3) + nl_state[a].l12 +
                nl_state[a].l3;
      int E = 2 * (nl_state[b].n12 + nl_state[b].n3) + nl_state[b].l12 +
              nl_state[b].l3;
      cout << "\t nl_p : " << a << "\t E_p = " << E_p << "\t nl : " << b
           << "\t E = " << E << endl;
    }
  } /* - */
  // ---------------------------------------------------------------------- //
  //
  // ------------------- Calculate PFKglbarS_vec ------------------- //
  cout << "  Calculate PFKglbarS_vec   " << endl;
  this->PFKglbarS_build();
  cout << "  Finished Calculate PFKglbarS_vec  " << endl;
  /* - */
}

void Wave_Fun::PFKglbarS_build() {

  int size_nlpnl = nlpnl_state.size();
  int size_nl = nl_state.size();

  // if (WaveCom.rank() == 0) {
  // cout << "size_nlpnl = " << size_nlpnl << "\t size_nl = " << size_nl
  //     << "\t on Rank : " << WaveCom.rank() << endl;
  //}
  // int size_lbar = lbar_state.size();
  PFKglbarS_vec.resize(4);
  int nunzero_num = 0;
  // cout << "Begin to init PFKglbarS" << endl;
  // Pgs_Bar pgs;
  for (int S_index = 0; S_index <= 3; S_index++) {
    // cout << "S_index : " << S_index << endl;
    // pgs.init(size_nl, 30);
    PFKglbarS_vec[S_index].resize(size_nlpnl);
    if (WaveCom.rank() == 0) {
      cout<<" S_index = "<<S_index<<"  size_nlpnl : "<<size_nlpnl<<endl;
    }
    for (int m = 0; m < size_nlpnl; m++) {
      // cout << "num : " << num << endl;
      if (WaveCom.rank() == 0) {
        // cout << "=== m = " << m;
      }
      int l_index = nlpnl_2_lpl[m];
      if (WaveCom.rank() == 0) {
        // cout << "l_index = " << l_index << endl;
      }
      int size_lbar = Llbar_Ic_Reduced[l_index].size();
      PFKglbarS_vec[S_index][m].resize(size_lbar);
      int num_p = nlpnl_state[m].nl_p;
      int num = nlpnl_state[m].nl;
      if (WaveCom.rank() == 0) {
        cout << "\t num_p : " << num_p << "\t "
             << "num : " << num << endl;
        cout << "\t size_lbar = " << size_lbar << endl;
      }
      //#pragma omp parallel for
      for (int lb_r = 0; lb_r < size_lbar; lb_r++) {

        int lb = Llbar_Ic_Reduced[l_index][lb_r];
        if (WaveCom.rank() == 0) {
          //  cout << "\t \t lb_r : " << lb_r << "\t lb : " << lb << endl;
        }
        int lbar = lbar_state[lb].lbar;
        int lbar12 = lbar_state[lb].lbar12;
        int lbar3 = lbar_state[lb].lbar3;
        int size_klamb = klamb_state[lbar].size();
        // cout<<"size_klamb : "<<size_klamb<<endl;
        PFKglbarS_vec[S_index][m][lb_r].resize(size_klamb);
#pragma omp parallel for schedule(dynamic, 1)
        for (int kl = 0; kl < size_klamb; kl++) {
          int A = klamb_state[lbar][kl].lamb12_p;
          int B = klamb_state[lbar][kl].lamb12;
          int C = klamb_state[lbar][kl].lamb3_p;
          int D = klamb_state[lbar][kl].lamb3;
          double T1 = PFKglbarS_vec_cal2(num_p, num, A, B, C, D, S_index, lbar,
                                         lbar12, lbar3);
          /*double T2=PFKglbarS_vec_cal(num_p, num, A, B, C, D, S_index,
          l); if(std::abs(T1 - T2 )>0.00000001){
              cout<<num_p<<"\t"<<num<<"\t"<<A<<"\t"<<B<<"\t"<<C<<"\t"<<D<<"\t";
              cout<<S_index<<"\t"<<l<<endl;
              exit(0);
          }*/
          PFKglbarS_vec[S_index][m][lb_r][kl] = T1;
          nunzero_num++;
        }
      }
    }
    //}
    // pgs.EndBar();
  }
  cout << "Total non-zero PFKglbarS_vec number = " << nunzero_num
       << "\t on rank " << WaveCom.rank() << endl;
}

double Wave_Fun::PFKglbarS_quick_t(int &num_p, int &num, int &A, int &B, int &C,
                                   int &D, int &S_index, int &lbar, int &lbar12,
                                   int &lbar3) {
  double res = 0.0;
  int lbar_index = this->lbar_vec[lbar][lbar12][lbar3];
  int klamb_index = this->klamb_vec[lbar][A][B][C][D];
  if (klamb_index < 0) {
    cout << "Wrong at PFKglbarS_quick: klamb_index <0 " << endl;
    cout << "A = " << A << "\t B = " << B << "\t C = " << C << "\t D = " << D
         << endl;
    exit(0);
  }
  int nlpnl = this->nlpnl_vec[num_p][num];
  if (nlpnl < 0 || nlpnl > this->nlpnl_2_lpl.size()) {
    // if (l_index < 0 || l_index > this->Llbar_Ic_Red_vec.size()) {
    std::cout << "nlpnl < 0 || nlpnl > this->nlpnl_2_lpl.size() on Rank : "
              << WaveCom.rank() << std::endl;
    cout << "num_p = " << num_p << "\t num = " << num << "\t nlpnl = " << nlpnl
         << endl;
    exit(0);
  }
  int l_index = this->nlpnl_2_lpl[nlpnl];
  if (l_index < 0) {
    // if (l_index < 0 || l_index > this->Llbar_Ic_Red_vec.size()) {
    std::cout << "Wrong at PFKglbarS_quick : l_index < 0 || l_index > "
                 "Llbar_Ic_Red_vec.size()"
              << std::endl;
    cout << "num_p = " << num_p << "\t num = " << num << "\t nlpnl = " << nlpnl
         << endl;
    exit(0);
  }

  if (lbar_index < 0 || lbar_index > Llbar_Ic_Red_vec[l_index].size()) {
    cout << "Wrong at PFKglbarS_quick : lbar_index < 0 || lbar_index > "
            "Llbar_Ic_Red_vec[l_index]"
         << endl;
    cout << "l_index = " << l_index << "\tLlbar_Ic_Red_vec[l_index].size() = "
         << Llbar_Ic_Red_vec[l_index].size() << endl;
    cout << "lbar_index = " << lbar_index << "\t :  lbar = " << lbar
         << "\t lbar12 = " << lbar12 << "\t lbar3 = " << lbar3 << endl;
    exit(0);
  }
  int lb_reduce = Llbar_Ic_Red_vec[l_index][lbar_index];
  if (nlpnl < 0) {
    cout << "Wrong at PFKglbarS_quick: nlpnl <0 " << endl;
    cout << "num_p = " << num_p << "\t num = " << num << "\t nlpnl = " << nlpnl
         << endl;
    exit(0);
  } /* - */
  res = PFKglbarS_vec[S_index][nlpnl][lb_reduce][klamb_index];
  return res;
}

double Wave_Fun::PFKglbarS_vec_cal2(int num_p, int num, int A, int B, int C,
                                    int D, int S_index, int lbar, int lbar12,
                                    int lbar3) {
  int n12_p = nl_state[num_p].n12;
  int n3_p = nl_state[num_p].n3;
  int l12_p = nl_state[num_p].l12;
  int l3_p = nl_state[num_p].l3;

  int n12 = nl_state[num].n12;
  int n3 = nl_state[num].n3;
  int l12 = nl_state[num].l12;
  int l3 = nl_state[num].l3;

  // int size_k123 = k123_state.size();
  // int k12, k12_p, k3, k3_p;
  // double part, part1;
  // int lbar = lbar_state[lbar_index].lbar;
  // int lbar12 = lbar_state[lbar_index].lbar12;
  // int lbar3 = lbar_state[lbar_index].lbar3;

  int lbar_index_g = lbar_vec_g[lbar][lbar12][lbar3];

  vector<double> PF(size_k123);
  vector<double> PF_p(size_k123);
  //#pragma omp parallel for
  for (int K = 0; K < size_k123; K++) {
    int k12 = k123_state[K].n_k12;
    int k3 = k123_state[K].n_k3;
    double F = F_reg_vec[k12][k3];
    PF_p[K] =
        F * Pwave_vecwx[A][n12_p][l12_p][k12] * Pwave_vecwx[C][n3_p][l3_p][k3];
    PF[K] = F * Pwave_vecwx[B][n12][l12][k12] * Pwave_vecwx[D][n3][l3][k3];
  }
  double part = 0.0;
  for (int K_p = 0; K_p < size_k123; K_p++) {
    double part1 = 0.0;
    for (int K = 0; K < size_k123; K++) {
      part1 += glbarS2_vec[S_index][lbar_index_g][K_p][K] * PF[K];
    }
    part += part1 * PF_p[K_p];
  }
  return part;
}

void Wave_Fun::evaluate() {
  // this function is to evaluate the memory and time cost for Wave
  // For PFKglbar0 & PFKglbar2
  int size_num = this->nl_state.size();
  int size_lamb1 = 3; // 0,1,2
  int size_lbar = this->lbar_max + 1;
  vector<double> Test(10);
  for (int i = 0; i < 10; i++) {
    Test[i] = 0.01;
  }
  double Test_t = 1.01;
  double mem_per_vec = sizeof(Test_t);
  cout << "size of a vector with 10 element : " << sizeof(Test) << endl;
  cout << " size_num = " << size_num << endl;
  cout << " size_lbar = " << size_lbar << endl;

  int size_PFKglbar0 = size_num * size_num;
  int size_PFKglbar2 = size_num * size_num * size_lamb1 * size_lbar;

  // For PFKglbarS_vec & glbarS2_vec
  int size_S = 4;
  int size_nlpnl = nlpnl_state.size();
  int size_lbar_state = lbar_state.size();
  int size_lbar_state_g = lbar_state_g.size();
  int size_k123 = k123_state.size();
  cout << " mem_per_vec = " << mem_per_vec << " bit" << endl;
  cout << " size_nlpnl = " << size_nlpnl << endl;
  cout << " size_lbar_state = " << size_lbar_state << endl;
  cout << " size_lbar_state_g = " << size_lbar_state_g << endl;

  int size_glbarS = size_S * size_lbar_state_g * size_k123 * size_k123;

  double size_PFKglbarS = 0;
  for (int l = 0; l < lbar_state.size(); l++) {
    int lbar = lbar_state[l].lbar;
    int size_klamb_state = klamb_state[lbar].size();
    size_PFKglbarS += size_klamb_state;
  }
  size_PFKglbarS = size_PFKglbarS * size_S * size_nlpnl;

  cout << "PFKglbar0_vec size : " << size_PFKglbar0 << endl;
  cout << "\t with Memory ~ " << size_PFKglbar0 * mem_per_vec / 1024.0 / 1024.0
       << "MB" << endl;
  cout << "PFKglbar2_vec size : " << size_PFKglbar2 << endl;
  cout << "\t with Memory ~ " << size_PFKglbar2 * mem_per_vec / 1024.0 / 1024.0
       << "MB" << endl;
  cout << "glbarS_vec size : " << size_glbarS << endl;
  cout << "\t with Memory ~ " << size_glbarS * mem_per_vec / 1024.0 / 1024.0
       << "MB" << endl;
  cout << "PFKglbarS_vec size : " << size_PFKglbarS << endl;
  cout << "\t with Memory ~ " << size_PFKglbarS * mem_per_vec / 1024.0 / 1024.0
       << "MB" << endl;

  double size_tot1 = size_PFKglbar0 + size_PFKglbar2 + size_PFKglbarS;
  cout << "PFKglbar0_vec + PFKglbar2_vec + PFKglbarS_vec size :" << size_tot1;
  cout << "\t with Memory ~ " << size_tot1 * mem_per_vec / 1024.0 / 1024.0
       << "MB"
       << "\t " << size_tot1 * mem_per_vec / 1024.0 / 1024.0 / 1024.0 << "Gb"
       << endl;
  /*
    size_PFKglbar0 = 0;
    size_PFKglbar2 = 0;
    size_glbarS = 0;
    size_PFKglbarS = 0;

    for (int i = 0; i < PFKglbar0_vec.size(); i++) {
      size_PFKglbar0 += PFKglbar0_vec[i].size();
    }
    for (int i = 0; i < PFKglbar2_vec.size(); i++) {
      for (int j = 0; j < PFKglbar2_vec[i].size(); j++) {
        for (int k = 0; k < PFKglbar2_vec[i][j].size(); k++) {
          size_PFKglbar2 += PFKglbar2_vec[i][j][k].size();
        }
      }
    }
    for (int i = 0; i < PFKglbarS_vec.size(); i++) {
      for (int j = 0; j < PFKglbarS_vec[i].size(); j++) {
        for (int k = 0; k < PFKglbarS_vec[i][j].size(); k++) {
          /*if (i == 0 && j == 0)
            cout << i << "\t " << j << "\t " << k << "\t "
                 << PFKglbarS_vec[i][j][k].size() << endl;
          size_PFKglbarS += PFKglbarS_vec[i][j][k].size();
        }
      }
    }
    for (int i = 0; i < glbarS2_vec.size(); i++) {
      for (int j = 0; j < glbarS2_vec[i].size(); j++) {
        for (int k = 0; k < glbarS2_vec[i][j].size(); k++) {
          size_glbarS += glbarS2_vec[i][j][k].size();
        }
      }
    }
    cout << "PFKglbar0_vec size : " << size_PFKglbar0 << endl;
    cout << "\t with Memory ~ " << size_PFKglbar0 * mem_per_vec / 1024.0 /
    1024.0
         << "MB" << endl;
    cout << "PFKglbar2_vec size : " << size_PFKglbar2 << endl;
    cout << "\t with Memory ~ " << size_PFKglbar2 * mem_per_vec / 1024.0 /
    1024.0
         << "MB" << endl;
    cout << "glbarS_vec size : " << size_glbarS << endl;
    cout << "\t with Memory ~ " << size_glbarS * mem_per_vec / 1024.0 / 1024.0
         << "MB" << endl;
    cout << "PFKglbarS_vec size : " << size_PFKglbarS << endl;
    cout << "\t with Memory ~ " << size_PFKglbarS * mem_per_vec / 1024.0 /
    1024.0
         << "MB" << endl;

    size_tot1 = size_PFKglbar0 + size_PFKglbar2 + size_PFKglbarS;
    cout << "PFKglbar0_vec + PFKglbar2_vec + PFKglbarS_vec size :" << size_tot1;
    cout << "\t with Memory ~ " << size_tot1 * mem_per_vec / 1024.0 / 1024.0
         << "MB" << endl;
        /* - */
  //
}
